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ABSTRACT 


The aim of the work is to investigate the aerodynamic characteristics such as lift coefficient, drag coefficient, pressure distribution over 

a surface of an airfoil of NACA-4312. A commercial software ANSYS Fluent was used for these numerical simulations to calculate 

the aerodynamic characteristics of 2-D NACA-4312 airfoil at different angles of attack (a) at fixed Reynolds number (Re), equal to 

5 x 10 * * * * 5 . These simulations were solved using two different turbulence models, one was the Standard k — 8 model with enhanced wall 
treatment and other was the SST k — co model. Numerical results demonstrate that both models can produce similar results with little 
deviations. It was observed that both lift and drag coefficient increase at higher angles of attack, however lift coefficient starts to reduce 
at a =13° which is known as stalling condition. Numerical results also show that flow separations start at rare edge when the angle of 
attack is higher than 13° due to the reduction of lift coefficient. 


Keywords: Airfoil; CFD; RANS; Lift Coefficient; Drag Coefficient. 
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1. Introduction 

Aerodynamics is defined as a study of the resulting effects of 
relative motion between air molecules and body surfaces [1]. 
Understanding the air flow behavior around a moving body is 
very important because we can control the air flow to our 
advantage in the case of aircraft, wind turbine, drone, and many 
more. Aerodynamics focuses on studying this phenomenon by 
applying basic physics laws such as Newton's laws and the 
Navier-Stokes equation. To study the aerodynamic profile of a 
particular object, it is necessary to define the object's shape, 
where the airfoil or the wing cross-section is used in the case of 
an airplane. The cross-section of the typical aircraft wing is an 
airfoil and is largely responsible for producing the 
forces that sustain the aircraft in flight [2]. Airfoil is known as 
the cross-sectional shape of a wing, blade (of a propeller, rotor, 
or turbine) which is placed in an airstream in order to generate 
useful aerodynamic forces. An airfoil's aerodynamic profile can 
be examined in two ways, one experimentally and other 
computationally. The aerodynamic characteristics of an airfoil 
are obtained in the experimental method by placing the airfoil (or 
wing model) in a controlled wind tunnel and recording the 
velocity and pressure distribution around the airfoil. The same 
experiment can be done in the computational approach using a 
CFD code. Both approaches are usually complementary to each 
other and can help us design a specific purpose for the airfoil [3]. 

In the early 1800s, it was discovered that a curved surface 
produces more lift than a similar size flat plate by Sir George 
Cayley [4]. The most effective way to do this is to use an airfoil. 
In the Langley two-dimensional low turbulence tunnels tests a 
considerable amount of airfoil data has been accumulated by 
National Advisory Committee for Aeronautics (NACA) in the 
United States. The development of NACA airfoils which are now 
in common use was started in 1929 with a systematic 
investigation of a family of airfoils in the Langley variable- 
density tunnel [5]. Airfoils of this family are designated by 
numbers having four digits. Experiments were performed to test 


different airfoils and to determine the lift and drag coefficient for 
different airfoils and to report their results. Jacobs et al. (1933) 
performed a wide variety of experiments comparing the 
geometric airfoil parameters and angle of attack and wind speed 
for different NACA airfoils and published the results of the lift 
and drag coefficient for 78 different NACA airfoils [6]. Since 
these experiments were carried out in a wind tunnel with a 
relatively high free-stream turbulence intensity of about 2%, the 
results are expected to vary from low free-stream turbulence 
conditions. NACA airfoils have been used extensively as 
aerodynamic test models. NACA-0012 and NACA-4412 
profiles are arguably the most studied airfoils. Ravi et al. [7] 
predicted numerically a transition model of an incompressible 
laminar to turbulent flow over NACA-4412 airfoil at Reynolds 
number of 3 x 10 6 . Eleni et al. [8] also carried out studies on 
the variation of lift and drag coefficients from flow around 
NACA-0012 airfoil at Reynolds number of 3 X 10 6 for different 
turbulence models. In addition to the results available for a 
NACA-0012 and NACA-4412 airfoil profile, a substantial 
experimental or numerical database is not available for NACA- 
4312 airfoil. This study is motivated by the need to complement 
the presently limited body of knowledge for NACA-4312 airfoil. 
Thus, the present investigation is focused on examining the effect 
of the Reynolds number and the angle of attack on the 
performance characteristics of a NACA-4312 airfoil and relating 
the performance characteristics. 

In this study, CFD approach is used to determine the 
aerodynamic profiles of the NACA-4312 airfoil. The CFD 
packages contain three main elements which are Pre-processor, 
Solver and Post-processor [9]. Recently CFD has been the 
method of choice in the aerospace, automotive and many 
industrial components. CFD is vastly used in the field of 
aerodynamics because it is cheaper than experimental process 
and also gives more accurate results. The transition from the 
laminar to the turbulent flow phase plays a very important role in 
simulating the flow over an airfoil. The fluid flow over the airfoil 
exerts a pressure force perpendicular to the upper and lower 
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surfaces along with shear force which is parallel to the surfaces. 
The resultant of these two forces is the area of interest. The 
normal component of the force is known as lift force and the 
force which is acting in the flow direction is known as drag [10]. 

It was observed from the previous literature review that there 
is a lack of understanding of the characteristics of NACA-4312 
airfoil. In this paper a series of numerical simulations were 
carried out to analyze the characteristics of NACA-4312 airfoil 
with cord length of 1000 mm and Reynolds number of 5 x 10 5 . 
ANSYS Fluent was used to solve the steady-state RANS 
(Reynolds Average Navier-Stokes) equation with two different 
turbulence models. 

2. Methodology 

2.1 Airfoil 

In this paper, NACA-4312 airfoil was selected for the 
simulations, the profile is shown in Fig. 1. The NACA-4312 
airfoil means that it has a maximum camber equal to 4% of the 
cord which is located at 30% of the cord from the leading edge 
with a maximum thickness of 12% of the chord. In this paper, 
cord length was considered 1000 mm. Airfoil was created by 
ANSYS design modeler by importing coordinate file of 
NACA-4312. 



0 0.2 0.4 0.6 0.8 1 

X(m) 


Fig. 1 Geometry of NACA-4312 airfoil. 

2.2 Computational Method 

These simulations were conducted by ANSYS Fluent. The 
problem was solved in steady-state with two turbulence 
models. These simulations were conducted at fixed Reynolds 
number (Re), equal to5xl0 5 , to show the effect at the 
transition region [11]. Air was assumed as working medium 
with a constant density (p ) and viscosity (ju), whereas p = 
1.225 kg/m 3 , and p = 1.7894 x 10 -5 kg/ms [12]. The steady- 
state RANS (Reynolds Average Navier-Stokes) equation was 
solved using the Green-Gauss cell based gradient option and 
pressure based solver was selected as the flow is 
incompressible. The RANS equations are time- 
averaged equations of motion for fluid flow, as bellows: 

— (pu i ) = 0 (1) 


—(P u iUj )= 


dp 8 

dx- dx - 

1 J 


CM; du • 
[jU(— + — 
dx - dx- 
J i 


2 8u, 
-<%—)] + 

3 « a, 


— {-pu'p'j) 


( 2 ) 


where p indicates density, u indicates velocity and p 
indicates dynamic viscosity of the fluid. The left side of this 
equation describes the fluid element's mean momentum change 
to the flow instability and convection by the mean flow. This 
adjustment is regulated by the mean body force, the isotropic 
stress due to the mean pressure field, the viscous stresses, and 
visible stress due to the fluctuating velocity field, commonly 
referred to as the Reynolds stress (—pu[u)). This nonlinear 
stress term requires additional modeling to solve the RANS 
equations and has resulted in many different turbulence models 
[13]. 

In this paper, the Standard k — £ model and the SST k — 
co model are used to predict the effects of turbulence in flow 
over the airfoil. 

2.2.1 Standard k-s Model 

K-epsilon (k — s) turbulence model is the most common 
model used in Computational Fluid Dynamics (CFD) for 
simulating flow characteristics for turbulent flow conditions. It 
is a two-equation model, which provides a general definition of 
turbulence via solving two transport equations (PDEs) one for 
the turbulence kinetic energy (k) and the other is its dissipation 
rate (£) [13]. The turbulent kinetic energy k and dissipation £ 
are obtained from the following transport equations (3) and (4): 

a 

- (pku -) = 

dx- 

8 u 8k (3) 


d d ju f ds 

- (pen- ) =- V(jU + ^- L ) -] 

dx- dx - cic- dx ■ 

1 J J 2 

+C 1 sk^ G k + C 3 s G b ) “ C 2s p ~k +S s 


(4) 


In these equations, 

Gk= Generation of turbulence kinetic energy due to mean 
velocity gradient 

YM=Fluctuation in compressible turbulence to overall 

dissipation rate 

C l£ , C 2£ , C 3e = Constants 

o kt o £ - Turbulent Prandtl number for k and £ 

S k ,S £ = User defined source terms 

k 2 

Here turbulent viscosity p t = pC M — ; where is C M a 
constant. 

2.2.2 SST k-co Model 

In CFD analysis SST (Shear Stress Transport) k — co is a 
widely used and two-equation eddy-viscosity turbulence 
model. It provides a general solution for the turbulence kinetic 
energy (k) and specific dissipation rate of eddy viscosity (co). 
This model gives more accurate results in numerical flow 
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analysis of airfoil than the standard k — co turbulence model as 
it includes transportation of turbulence sheer stress. The proper 
transport behavior can be obtained to the formulation of eddy- 
viscosity [14] .The turbulent kinetic energy k and dissipation 
rate of eddy viscosity co are obtained from the following 
transport equations (5) and (6): 

± tt , t ^± Kll ^)^G k -Y k+ S k (5) 

^ j j 

d d dco 

— (pcoitj ) = — [(// + ) — ] + Gfo - + Sfo (6) 

ox- ox - ox - 

1 J J 


In these equations, 

Gk= Generation of turbulence kinetic energy due to mean 
velocity gradient 

= Generation of co 
Y k ,Yu = The dissipation of k and co 
a k ,c 7 o) = Turbulent Prandtl number for k and co 
Sk’Scj - User defined source terms 


Here turbulent viscosity 


pk 

00 


r 1 SF? n ! 
max ———\ 
l a*’a 1(]d x 


where S is the 


strain rate magnitude, F 2 is a constant, and 
oc=Angle of attack. 


2.3 Boundary Condition with domain 

For these simulations, a computation domain was created 
around the NACA-4312 airfoil, shown in Fig. 1. Chord length 
of NACA-4312 airfoil was assumed to be equal to 1000 mm. 
To minimize the boundary effect, the domain was extended 
12.5C (C indicated chord length) along the upstream and 20C 
along downstream from the trailing edge, like Fig. 2 [15]. In 
the domain, airfoil surface was assumed no-slip boundary 
condition, BAFED was assumed constant velocity inlet and 
pressure outlet at BCD. For this fixed Re , the inlet velocity u, 
was assumed to 7.5 m/s. In these simulation angle of attack (a) 
was changed by changing flow direction instead of rotating 
airfoil which produced the same effect on the airfoil. The 
velocity components along v and y direction were calculated by 
u cos a and u sin a, respectively, for different angle of attack 
(a). 


2.4 Mesh generation and Wall Treatment 

The C-type structural mesh was created for better 
convergence and control of the wall function, shown in Fig. 3 
and Fig. 4. Fig. 3 depicts the structural meshing with 
quadrilateral elements which was done by ANSYS Meshing 
whereas Fig. 4 illustrates the meshing quality around the airfoil 
profile. 

The application of wall functions near the wall region may 
significantly reduce the processing and storage requirements 
while producing an acceptable degree of solution. To ensure 
sufficient boundary layer modeling the inflation was set to 1.15. 

The effectiveness of k — s model is beyond y + = 30. The 
non-dimensional wall parameter is defined as: 


A Velocity Inlet B 




Fig. 3 Mesh among the whole domain. 



Fig. 4 Mesh around the NACA-4312 airfoil. 
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+ yj{rcol p) 

y =y* - 

A 

Here y is the distance from the wall to the centroid of the 
first fluid cell. It was observed for all simulations that the value 
of y + was higher than 30 and lower than 60. 

2.5 Mesh Independence Test 

A series of simulations were conducted for mesh 
independency tests. A different number of meshes were created 
by diving circular and rectangular sections with different 
numbers. Fig. 5 depicts the variation of lift coefficient for 
different numbers of elements of the mesh for an angle of attack 
5°. It is observed from figure that meshes with higher than 
105000 number of elements can produce accurate results with 
minimum deviation. As a result mesh with element number 
105000 were selected for further simulation. 



Fig. 5 Variation of coefficient of lift with no. of element. 

3. Results and Discussion 

Fig. 6 depicts the variation of lift coefficient for different 
angles of attack for two different turbulence models. It is 
observed that the lift coefficients are almost equal for both 
models up to a = 6°, however, the SST k — co model predicts 
higher lift coefficient compared to the standard k — s model 
when a > 6°. Fig. 6 also shows that the lift coefficient increases 
with angle of attack, for both models, up to 12°. Later on, it 
decreases with a due to the flow separation at the trailing edge, 
which indicated stalling occurs at the range of 12° - 13°. From 
0° angle of attack to 12° angle of attack the lift curve is almost 
linear. 

Fig. 7 shows the effect of the angle of attack (a) on the 
drag coefficient for two different turbulence method. It is 
observed that the drag coefficient estimates higher values with 
the standard k — £ model compared to the S ST k — co because 
of the transportation of the shear viscosity. For both models, the 
drag coefficient increase with a , however, the increment rate is 
higher at higher a due to the flow separation at trailing edge. 
The drag coefficient still increases after the stalling angle 
however lift coefficient decreases. 



Fig. 6 Lift coefficient vs angle of attack. 



Fig. 7 Drag coefficient vs angle of attack. 



Fig. 8 Pressure coefficient on the airfoil surface at 10° angle 
of attack. 
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Fig. 9 Contours of pressure with the Standard k-s model for different angles of attack of NACA-4312 airfoil. 
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Fig. 10 Contours of pressure with the SST k-co model for different angles of attack of NACA-4312 airfoil. 
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Fig. 11 Velocity contours with the Standard k-e model 
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Fig. 12 Velocity contours with the SST k-co model 
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Fig. 8 demonstrates pressure distribution on the airfoil 
surface for both models at an angle of attack a - 10°. The upper 
and lower portion of the graph represent pressure coefficient 
distribution on the lower surface and upper surface of the 
airfoil, respectively. It is observed that the upper surface have 
lower pressure and lower surface have higher pressure which 
produces lift. Maximum lower and higher pressure observed 
near to the leading of the lower and upper surface, respectively, 
which starts to reduce along the trailing edge. 

3.1 Pressure Contours 

Fig. 9 and Fig. 10 presents pressure contours for the 
Standard k — £ turbulence model and the SST k — co 
turbulence model, respectively, at different angle of attack. In 
these two figures, the effect of angles of attack is shown by 
changing the flow directions instead of rotating airfoil, which 
gives the same effect on airfoil. Similar pressure distribution 
pattern is observed for the Standard k — £ and SST k — a) 
turbulence model, present in Fig. 9 and Fig. 10, respectively at 
lower angle of attack, whereas slightly varies in magnitude at 
higher angle of attack though the pattern of pressure 
distribution is similar. It is clearly observed from both figures 
that the upper surface has lower pressure and the lower surface 
has higher pressure. It is also observed from both figures that 
negative pressure region on the upper surface which is situated 
almost around the whole surface which starts to decrease with 
the increase of angle of attack, however pressure on the lower 
surfaces increases with angle of attack. Though pressure 
difference between lower and upper surface increase until flow 
separates on the upper surface at a certain angle of attack. The 
lift coefficient increases in both turbulence model with the 
increase of angle of attack and it reaches the stalling position 
near 13° angle of attack. Between the two turbulence models, 
it is observed that the pressure difference in the SST k — co 
model is slightly higher than the Standard k — s model which 
is clearly observed in Fig. 8 also. 

3.2 Velocity Contours 

Fig. 11 and Fig. 12 depict velocity contours for the 
Standard k — e turbulence model and SST k — co turbulence 
model, respectively, at different angles of attack. In these two 
figures, the effect of angle of attack is shown by changing the 
flow directions instead of rotating airfoil, which gives the same 
effect on airfoil. Similar velocity distribution pattern is 
observed for the Standard k — £ and SST k — co turbulence 
model, present in Fig. 11 and Fig. 12, respectively at lower 
angle of attack, whereas slightly varies in magnitude at higher 
angle of attack though patter of velocity distribution is similar. 
It is clearly observed from both figures that velocity above the 
upper is higher than the velocity below the lower surface. This 
happens due to the Bournollie’s principle. It is also observed 
from both figures that at lower angle of attack flow is attached 
to upper surface and lower surface as a result there is no zone 
of zero velocity which indicates flow separation happens. 
When the angle attack reaches to the value of 10° flow starts to 


separate, see Fig. 11, at the trailing edges which starts to move 
towards the leading edges at higher angles of attack. A similar 
pattern is also observed for the SST k — co model, see Fig. 12, 
however, flow separation is identified comparative lower angle 
of at 8° compare to the Standard k — £ model. Due to this flow 
separation, lift coefficient starts to decrease after a certain angle 
of attack. 

4. Conclusion 

In this paper aerodynamic characteristics of NACA-4312 
airfoil have been studied using ANSYS Fluent by the Standard 
k — £ model and SST k — a) model turbulence model. At lower 
angles of attack both methods produce similar results however 
at higher angles of attack results vary in magnitude slightly 
though nature is similar. It is observed that the lift and drag 
coefficient increase with angle of attack though the lift 
coefficient shows a decreasing trend beyond the angle of attack 
13°. 
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ABSTRACT 

Novel luffa fiber reinforced epoxy composites are prepared and their mechanical properties are investigated before and after chemical 
treatment. The unique natural knitting structure of luffa provides an excellent reinforcement to the epoxy matrix. Knowing that the 
fiber-matrix bond gets stronger and imparts more strength to the composite when chemical treatment is done on fibers, composites are 
manufactured by untreated and treated luffa fiber using epoxy as a matrix. Luffa fiber is treated using benzoyl chloride and NaOH. 
Tensile and flexural tests are conducted on composites to investigate the effect of chemical treatment. Test results have shown that the 
chemical treatment on fibers improved the tensile strength, tensile modulus, flexural strength and flexural modulus by 27.21%, 49.37%, 
41.84% and 6.44% respectively. Tensile modulus of luffa fiber composite is found to be higher compared with commonly used natural 
fiber composites. The experimental investigation suggests that, chemically treated luffa fiber reinforced epoxy composites could be a 
potential lightweight material in various applications. 


Keywords: Luffa Fiber; Benzoyl Chloride Treatment; Mechanical Properties; Natural Fiber Composite 



This work is licensed under a Creative Commons Attribution-NonCommercial 4.0 International 


1. Introduction 

In the recent time, composite materials have gained huge 
attention to the researchers due to its superior mechanical 
properties like strength to weight proportion, crack durability, 
better physical, chemical and electrical properties etc. [l]-[3]. In 
composite materials, fibers are used as reinforcement while the 
matrix conforms the shape of the material which also transfers 
the load between fibers. Fibers could be synthetic (glass, carbon, 
nylon etc.) or natural (jute, bamboo, luffa etc.). Although, 
synthetic fibers have better strength, natural fiber are 
comparatively cheaper, abundant in nature and renewable 
sources of fiber [4]. Despite the fact that, the strength of the 
natural fiber is lower than the synthetic fiber, natural fibers are 
used in automobile industries and many other products like ball 
point pen, fishing hook, tennis racket, bicycle etc. [5]-[6]. 

Previous researchers have studied mechanical properties of 
different natural fiber reinforced composites with or without 
fiber surface treatment. Akil et al. studied the water absorption 
of jute fiber reinforced polyester composites [7]. A.V.R. Prasad 
and K.M. Rao studied mechanical properties of jo war, sisal and 
bamboo fiber reinforced composites [8]. A lot of research has 
been carried out on bamboo fiber reinforced composites by 
different researchers [9]-[13]. Moisture absorption 
characteristics of sisal and roselle fiber reinforced hybrid 
composites were studied by Athijayamani et al. [14]. Dhakal et 
al. [15] studied the effect of water absorption on the mechanical 
properties of hemp fiber reinforced composites. Palm and coir 
fiber reinforced composites were studied by Haque et al. [16]. 

Luffa cylindrica is a natural fiber of less weight, low cost 
and useful mechanical property but also has an uncommon 


woven form which is exceptional in other natural fibers. 
Although a lot of research about mechanical properties of natural 
fiber reinforced composites were found in the literature. 
However, mechanical properties of luffa fiber reinforced 
composites are scarce in literature. Different mechanical 
properties of luffa fiber reinforced composites were studied by 
N. Mohanta and S.K. Acharya [17]-[18]. In this paper, the tensile 
and flexural properties of natural Luffa fiber are reported. The 
effect of NaOH and benzoyl chloride treatment on fiber surface 
to enhance the mechanical properties of Luffa fiber reinforced 
polymer composite with four layers of fiber interface is analyzed. 

2. Experimental Details 

2.1 Materials 

Luffa fiber is used as reinforcement and epoxy (LY 
556/HY 951 hardener) is used as matrix material. The Luffa 
cylindrica is collected from local agricultural market of 
Bangladesh and later, this luffa cylindrica is used to extract 
luffa fiber as described in section 2.2. The hardener is mixed 
with resin at a ratio of 1:10 as suggested by the manufacturer. 
NaOH and Benzoyl chloride are used to treat the fiber surface. 

2.2 Fiber Extraction 

Luffa fiber is collected out of ripe Luffa (Fig. 1(a)). When 
Luffa is ripe the xylem fiber network is opened consisting of 
lignin and cellulose. The natural lignin-cellulose bond is 
weaker. To form a stronger bond lignin needs to separate out 
which acts as matrix here. The fiber is separated from lignin- 
cellulose bond by fermenting in fermented water for 15 days 
(Fig. 1(b)). The extracted fiber is then dried in atmospheric air 
(Fig. 1(c)). 
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(a) (b) (c) 

Fig. 1 (a) Raw luffa, (b) Fermenting in water and (c) Drying in atmospheric air. 



(a) (b) 

Fig. 2 (a) Fiber mixing with NaOH solution and (b) After lh the NaOH and fiber mixture 





(a) (b) 

Fig. 3 (a) After 15 minutes continuous stirring the benzoyl chloride and fiber mix and (b) Drying in atmospheric air. 


2.3 Chemical Treatment 

To enact the hydroxyl bond of the cellulose and lignin, the 
strands are treated with 10% NaOH for 1 hour as shown in Fig. 
2. After that, the strands are kept in benzoyl chloride for 15 
minutes with continuous stirring (Fig. 3). Benzoyl chloride 
removes the cellulose and lignin bond by creating cellulose and 
benzoyl chloride bond. The strands are then soaked in ethanol 
which weakens the benzoyl chloride-cellulose bond. It is 
necessary to weaken this benzoyl chloride-cellulose bond 


because it reduces the fiber matrix interfacial strength. The 
strands are then washed properly using distilled water to 
remove any additional ethanol. Finally, the strands are dried in 
atmospheric air followed by oven drying at 70 °C for 6 hours 
[17]. 

2.4 Manufacturing Process 

Hand lay-up method is used which is a very effective way 
to manufacture composite [19]-[20]. Rectangular aluminum 
plate with a size of (300 mm x 300 mm) is used as mold plate. 
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Fig. 4 shows the schematic diagram of the lay-up. At first, the 
mold release is applied to the mold plate which facilitate the 
easy removal of composite after curing. The resin hardener 
mixer is then applied to the mold plate with a soft roller and 
brush. The first layer of fiber is then placed on the mold plate 
and applied resin hardener mixer with the brush and roller. The 
roller ensures uniform distribution of resin hardener mixer on 
the fiber layer and removes any air bubbles that formed. The 
next layer of fiber is then placed on top of the first layer and the 
procedure repeated. Top mold plate is placed on top of the top 
layer of fiber and a pressure of 100 MPa is applied which 
ensures the uniform thickness of the manufactured composite. 
At this stage, the composite is cured at room temperature for 24 
hours. In this way, composites are manufactured by using both 
treated and untreated fiber. Untreated fiber is used as a baseline 
for the comparison of mechanical properties. 



Fig. 4 Schematic diagram of lay-up. 


2.5 Tensile Test 

Tensile test on both type of composite (treated and 
untreated) are performed. Total 10 samples are tested among 5 
are from untreated fiber and the remaining are from treated 
fiber composite. The specimens are cut from the composite 
plate manufactured following the standard ASTM D3039 [21] 
and the specimen size is 147 mm x 20 mm x 4.5 mm as shown 
in Fig. 5. During the test, the specimens are placed inside the 
hydraulic grip of universal testing machine at a loading rate of 
2 mm/min. The load-displacement response are recorded and 
later these data are used to calculate stress and strain along with 
specimen width, thickness and a gauge length of 25 mm. 


elasticity, 0.001 to 0.003 strain range is selected as suggested 
by ASTM D7264. 


133 mm 



Fig. 6 Flexural test specimen. 

3. Results and Discussion 

3.1 Tensile Test Results 

As mentioned in section 2.5, tensile test is performed on 
luffa fiber reinforced epoxy composite for both untreated and 
treated fiber. Fig. 7(a) shows the stress-strain response for 
untreated fiber while Fig. 7(b) shows the same for the treated 
fiber. By comparing Fig. 7(a) and Fig. 7(b) it is observed that 
the maximum stress for treated fiber is higher compared to the 
untreated fiber which indicates that the strength of the 
composite manufactured from the treated fiber increases with 
fiber treatment. Fig. 8(a) confirms this statement, which is a 
bar chart showing the strength of untreated and treated fiber 
composite. The average tensile strength of untreated and treated 
fiber are 11.43 MPa and 14.54 MPa respectively, which 
suggests a 27.21% increase in the tensile strength as shown in 
Table 1. 



Strain (mm/mm) 

(a) 



2.6 Flexural Test 

Three-point bending test is performed on the manufactured 
specimen to evaluate the flexural properties. Five specimen 
from the untreated and five specimen from the treated fiber are 
tested following the standard ASTM D7264 [22] and the 
specimen dimension is 133 mm x 28 mm x 4.5 mm as shown 
in Fig. 6. Three-point bending tests are carried out on the same 
universal testing machine with a loading rate of 1 mm/min. A 
constant span length to thickness of 32:1 is maintained during 
the test and for calculating the flexural chord modulus of 



Fig. 7 Stress-strain diagram for the tensile test of (a) bamboo 
strip and (b) hybrid composite 
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Table 1 Tensile test results with standard deviation. 


Condition 

Tensile Strength (MPa) 

% Increase 

Tensile Modulus (GPa) 

% Increase 

Untreated 

11.43 + 0.81 

- 

2.37 + 0.41 

- 

Treated 

14.54+1.43 

27.21 

3.54 + 0.57 

49.37 



Untreated fiber 

(a) 


i Treated fiber 



Untreated fiber 

(b) 


i Treated fiber 


Fig. 


8 Effect of chemical treatment on (a) Tensile strength and 
(b) Tensile modulus. 


It is also observed from Fig. 7 that the slope of the stress- 
strain curve for treated fiber composite is higher compared to 
untreated fiber. This suggests that the tensile modulus increases 
with the treatment of fiber and is represented by Fig. 8(b). 
Table 1 shows a comparative analysis of the tensile modulus 
for untreated and treated fiber. The average tensile modulus for 
untreated fiber is 2.37 GPa and for treated fiber is 3.54 GPa. 
So, the treatment of fiber increases the tensile modulus by 
49.37%. 


is observed from this figure that the maximum strength and 
slope of the curve of treated fiber (Fig. 9(b)) is higher compared 
to the untreated (Fig. 9(a)) one. This is represented in Fig. 10 
and Table 2. Table 2 shows that the flexural strength and 
modulus of treated fiber is increased by 41.84% and 6.44% 
respectively compared to untreated fiber. 
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3.2 Flexural Test Results Fig. 9 Stress-strain diagram of flexural test for (a) Untreated 

fibers and (b) Treated fibers. 

Fig. 9 shows stress-strain response from the three point 
bending test for both untreated and treated fiber composites. It 


Table 2 Flexural test results with standard deviation. 


Condition 

Flexural Strength (MPa) 

% Increase 

Flexural Modulus (GPa) 

% Increase 

Untreated 

26.29 + 3.36 

- 

2.02 + 0.27 

- 

Treated 

37.29 + 7.09 

41.84 

2.15 + 0.61 

6.44 
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Untreated fiber ■ Treated fiber 


(a) 



Fig. 10 Effect of chemical treatment on (a) Flexural strength 
and (b) Flexural modulus. 


4. Conclusion 

The mechanical properties of luffa fiber reinforced epoxy 
composite are investigated in this paper. Luffa fiber is treated 
with NaOH and benzoyl chloride and the mechanical properties 
are compared for untreated and treated fiber composites. It is 
found that the tensile strength and modulus for treated fiber are 
increased by 27.21 % and 49.37 % respectively. Flexural 
strength and modulus is increased by 41.84 % and 6.44 % 
respectively. It is found that the overall mechanical properties 
of the composites improved because of chemical treatment. 
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ABSTRACT 


It is often a challenge to achieve uniform flow in turbulent swirl flow and to predict the flow within the nozzle as measurement 
diagnostics face difficulty to capture both mean flow and turbulence. The purpose of this study is to numerically investigate the near 
wall flow characteristics and turbulent behavior for the effect of different tangential inlet numbers of an incompressible turbulent 
swirl air jet. In this regard, axial-plus-tangential flow based swirling nozzle is considered for the simulation using finite volume 
method, where turbulence is approximated by the Shear Stress Transport (SST) k-co model. The results show that axial and tangential 
velocity at the wall vicinity response the most. Moreover, the turbulent flow characteristic for no swirl flow is nearly uniform, but 
for swirl flow it fluctuates abruptly near the inlet section where the swirl has introduced. The skin friction coefficient for 2TP is the 
maximum for swirl flow and for no swirl condition the skin friction coefficient is nearly uniform. Due to the swirl introduction the 
pressure drop characteristics near the nozzle center response quickly and near the wall vicinity this property changes slowly. The 
magnitude of swirl decay fluctuates before the nozzle converging section however after the nozzle converging section the swirl 
decay is nearly constant. The local swirl near the inlet is highly unpredictable although after the nozzle converging section the local 
swirl profile is nearly similar for 2TP, 3TP and 4TP. 


Keywords: Turbulent; Nozzle; Swirl; Aerodynamic; CFD; Pressure. 
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1. Introduction 

Swirl flows are a classical form of fluid flow and typically 
encountered in many engineering applications, such as gas 
turbine combustor and cyclonic separator, as well as in nature, 
such as tomedo. In swirling flows generally two types of vortex 
emerge, namely solid body rotation and free vortex flow. 
However, in actual practice, a mixed of these two typically 
appears whereby azimuthal velocity increases with radius in the 
core and decreases afterwards. The academic and industrial 
research on swirling jets was found to be challenging due to their 
unique characteristics, such as flow reversal, vortex bubble and 
anisotropic turbulences in both free and impinging jets. The 
detailed flow physics, instability and above features of swirling 
flows are explicitly available in the literature [1]. The axial and 
radial pressure gradient in swirling flows affects downstream 
flow development and wall-bounded flows. Swirl flows are 
generated in different applications in many ways, such as rotating 
pipe, rotating vanes and twisted tapes inside stationary pipe, 
axial-plus-tangential entry [2]-[7]. The wide variations of 
generations can be divided into two broad groups: geometrical 
and aerodynamical. In geometrically generated swirl flows, for 
example rotating vanes or twisted tapes inside pipes, the flow is 
disturbed by such geometries or obstacles, which exacerbates the 
generality of the flow and difficult to draw a summative 
observation. On the other hand, aerodynamically generated swirl 
flows, such as rotating pipe or axial-plus-tangential entry flow, 
have better control of the flow and common flow physics may 
establish in different researches. Despite geometrically generated 
swirl flows are widely available in the literature and still 
progressing, fundamentals of aerodynamic swirl flows are 
currently being investigated in detail because of its prevalence to 
use recently in number heat treatment applications. 


Yajnik and Subbaiah [8] experimentally studied the effects 
of swirl on internal turbulent flows by conducting experiments 
on turbulent pipe flow. They found that the mean velocity 
profiles close to and away from the wall admit similarity 
representations at sufficiently large Reynolds numbers and 
extended velocity defect law is sensitive to swirl as the wall law 
is not sensibly dependent on swirl. Chang and Dhir [9] found a 
flow reversal region in axial velocity profile in the central portion 
of an acrylic tube and an axial velocity increase near the wall. 
Kito and Kato [10] studied the near wall velocity distribution of 
turbulent swirling flow in circular pipes and concluded that the 
flow becomes three-dimensional after transitional swirl intensity. 
Later, Kitoh [11] investigated the turbulent behavior of free- 
vortex-type swirling flow through a long straight circular pipe 
and reported that the swirling component decays downstream as 
a result of wall friction. Buschmann et al. [12] explored the wall 
behavior, the location of the peak Reynolds shear stresses and the 
three normal stresses of turbulent channel/pipe flows, and stated 
that no scaling works equally for all parameters. Ahmed et al. 
[13] conducted measurements of mean velocity and turbulence 
of swirling flows using dual wire CTA probe both in the core and 
near-wall regions. However, the measurements were confined to 
the immediately above the nozzle exit plane. Effect of viscosity 
and surface tension of fluids and associated instabilities in 
annular flow were also analyzed in some studies [14]-[15]. 

Since both DNS and LES typically require sufficiently 
smaller sized mesh or large mesh quantity to resolve small-scale 
turbulence, Reynolds Averaged Navier-Stokes (RANS) were 
widely used in swirling flows by many investigators. Although it 
is believed that among various RANS models the Reynolds 
Stress Model (RSM) would perform better because of its ability 
to capture anisotropy of turbulent shear stress, this hypothesis 
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was not found to be universal in the past [16]-[17]. In fact, some 
studies revealed that two-equation models (based on Boussinesq 
approximation) can be applied with acceptable accuracy in 
comparison with experimental data for moderate swirl flows - 
[18]-[19]. Nouri-Borujerdi and Kebriaee [20] simulated the 
turbulent boundary layer of an incompressible viscous swirling 
flow through a conical chamber using finite volume method and 
reported that the boundary layer thickness is dependent on the 
velocity ratios, Reynolds number and nozzle angle. Najafi et al. 
[21] investigated turbulent swirling decay in a vertical straight 
fixed pipe where swirl in induced by rotating honeycombs. They 
showed RSM with two-layer zone model for different near wall 
approaches are fairly well to predict the swirling flow but fail to 
predict the pressure distribution along the pipe wall. Islam et al. 
simulated an aerodynamically generated swirl flow using SST k- 
co model in the nozzle exit plane. The result showed centerline 
velocity [22] decay for introducing low levels of swirl into the 
impinging jet and a significant reduction in turbulent kinetic 
energy at the wall region. More recently, researchers [23]-[27] 
performed simulations of aerodynamic swirl flow using SST k- 
cd model and validated their results with the literature. 



(a) 3D Model (fo )Different Tangential 

Ports 


The above discussion revealed that although a significant 
amount of research available on swirling flows, they are mostly 
confined to gas turbine combustors, cyclonic separator, or 
geometry induced swirl flows. Relevant aerodynamically 
swirling flows are either classical in nature (experimental) with 
outdated experimental facility to resolve near-wall and 
turbulence behaviors inside the nozzle or tested robustness of 
various numerical schemes against experimental data with 
limited focus on flow behaviors. Thus, a detailed study for the 
effect of geometric parameters on near-wall and turbulence 
behaviors for an aerodynamic swirl nozzle appears to be limited 
in the literature. As such, the current numerical study will bridge 
this gap by investigating non-swirling and swirling flows from 
an aerodynamic nozzle for the same initial and flow conditions. 
The paper will examine the effect of number of tangential ports 
on mean and turbulence flow development along the length of 
the nozzle. 

2. Methodology 

2.1 Problem formulation 

An aerodynamically generated swirl nozzle, which is 
capable of seamless transition from non-swirling to highly 
swirling jets, is considered in this study. The nozzle is axial- 
plus-tangential entry type and consists of three tangential 
around the nozzle periphery and an axial port at the bottom of 
the nozzle. Detailed dimension of this nozzle is available in 
[4],[28], hence is not repeated here for brevity. The uniform 
flow was ensured by the flow settling chamber with 
honeycomb. The aerodynamic swirl flow was generated when 
both the uniform axial and tangential flows from different 
number of circumferentially oriented and inclined ports mix 
together farther downstream. In this case, three variants of 
tangential ports: 2-Tangential ports (2TP), 3-Tangential ports 
(3TP) and 4-Tangential ports (4TP) are considered. The exit of 
the nozzle has diameter (. D ) 40 mm, with a total length of 577 
mm. A three-dimensional view with the relative orientation of 
tangential ports is shown in Fig. 1. 


Fig. 1 (a) CAD view of the aerodynamic swirl nozzle (sliced 
to show the internal cavity), and (b): different number of 
tangential ports to impart azimuthal component. 


2.2 Governing Equation and Boundary Conditions 


The governing equations to solve the incompressible, 
steady-state flow characteristics within the aerodynamic swirl 
nozzle are the conservation of mass and momentum as follows: 


v-y = 0 


a) 


pV-VV = -VP + /iVT 


( 2 ) 


Here, ^ is the velocity vector, and p and p are the density 
and dynamic viscosity of the fluid, respectively. Since the 
current problem is turbulent in nature, RANS approach is 
applied to solve the mean and turbulence quantities. In RANS 
approximation, each variable is composed of time-averaged 
part (steady) and turbulence part, as shown below: 


$-(/)+(/)' ^ 

where ^ is a variable used in equations (1) and (2). Upon 
implementing the equation (3) into the governing equations and 
using equation (1) and setting time-average of turbulence 
equals to zero, the resulting RANS equations emerge. The 
RANS equations are similar to the governing equations, except 
an additional term ( pu[Uj ) in equation (2). This additional term 
is known as turbulent shear stress, which governs the 
turbulence characteristics. The shear stress components are 
determined via mean velocity gradients by the Boussinesq 
hypothesis: 


u i u 'j =- k $ij - Mt 


8w dw ■ 


dx, 

v J 


- + - 


dx; 


(4) 
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Where, jut is the turbulent viscosity and is a function of k and co, 
which are determined via an appropriate turbulence model. In 
this research, the SST k-co model, proposed by Menter and 
Egorov [29], is chosen to model turbulence transport quantities: 
turbulence kinetic energy, k , and the specific dissipation rate, 
co. SST k -co model was found to be one of the best performing 
RANS models in swirling and wall-bounded flows [16], [30]- 
[31]. 

The boundary conditions considered in this problem are 
mass flow inlet at the axial and tangential ports, pressure outlet 
at the nozzle exit plane and wall at the nozzle wall. The inlet 
conditions are adopted from the study [4] for realistic 
predictions. For no-swirl flow, axial flows (only) are provided 
from the bottom of the nozzle, with no flows from tangential 
ports. Conversely, for the swirling flows, three tangential ports 
supply the same amount of mass flow, with no axial flow from 
the bottom of the nozzle. This ensures the same Reynolds 
number, where the average velocity in Reynolds number is 
determined by volume flow rate divided by nozzle exit area. 
The pressure outlet condition is applied at the nozzle exit with 
atmospheric pressure whereby turbulence is specified by 2% 
intensity and hydraulic diameter. Finally, no-slip condition at 
the nozzle wall with ambient temperature is used. 

2.3 Numerical Settings and Model Validation 

The above governing equations (1-4) and two transport 
equations (SST k-co) are solved by finite volume based 
commercial software ANSYS Fluent vl7. The pressure-based 
coupled algorithm is used to simultaneously solve the 
governing equations. For the pressure discretization PRESTO 
(PREssure STaggering Option) is applied, and the second-order 
central-differenced for diffusion terms and second-order 
upwind scheme for convective terms of the transport equations 
are used. Hexahedral mapped mesh type was used in multi¬ 
zone meshing with an element size 0.0025m. A typical mesh 
independence test is shown in Fig. 2 for 3TP. In this case, four 
different grid elements are tested, namely, 974001, 1310326, 
1647089 and 1892870 elements. The results are found to be 
invariant, except the 974001 elements. As such, a mesh 
containing 1647k element is chosen in this study. Another 
testing is also done for first layer height on y + values and wall 
shear stress profiles, and is shown in Fig. 3. It appears the first 
layer cell height has little influence on those wall 
characteristics. As such, a first layer height of 0.05 mm is 
chosen in this study. The figure also shows that above settings 
ensured y + values less than 0.2 for the whole domain. The 
solution is assumed to be converged when the residuals of the 
flow parameters are less than 10" 5 . Mass conservation is also 
checked between inlet ports and the nozzle exit for converged 
solutions and the difference between inlets and outlet is found 
close to zero. Cartesian coordinates are given by (v, y, z) with 
corresponding velocity components ( U x , U y , U z ). For 
convenience, results of this 3D simulation are presented by 
polar coordinates (v, r, 0) with x-axis coinciding for each 
system. In this regard, the corresponding axial (U) and 
tangential (W) velocity components are derived from the 
Cartesian coordinate data using the axes transformation rule. 



Fig. 2 Sensitivity of mesh elements for the case 3TP. 



(a) 



(b) 


Fig. 3 Effect of first layer cell height on y+ and normalized 
wall shear stress. 


sum of flowrates through all axial and tangential inlets. 
Mathematically, 

Q = Q 

Qa + Qt Or (5) 

The Reynolds number is defined as, 


2.4 Parameter Definition 


Av 


( 6 ) 


The flow ratio Q r is the ratio of total mass flow rate 
through tangential inlets to the total flow rates in the nozzle, i.e. 
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Where, A is the nozzle area at the exit plane of the nozzle 
and v is the dynamic viscosity of fluid. Two different swirl 
number definitions are common as the ratio of axial flux of 
tangential momentum to the axial flux of axial momentum: 


s = - 


2 jcp^ r 2 UWdr 


InRp^Au 2 -0.5W 2 \dr 


JO 


f r~UWdr 
o' _ *KJ_ 

rR 2 

In rU dr 


R 


Finally, the local swirl number is defined as, 


(7) 


( 8 ) 



* 


s 



2.5 Validation 


(9) 


The simulation data of the swirl nozzle is first tested with 
different turbulence models and validated by comparing with 
experimental data [4] at the nozzle outlet plane with three 
tangential ports for the case Q r = 0 and 1 (Fig. 4). It appears 
that SST k-co model predicts the flow behavior well for both 
non-swirling and swirling cases. As such, SST k-co model will 
be used for the data presented in the ensuing results and 
discussion section. Fig. 4 also shows a good agreement 
between the numerical prediction and the experimental data for 
Q r = 0. A slight deviation is observed for Q r = 1, but 
importantly, the numerical data predicts well the profile 
behaviors and peak locations. The deviation is attributed to the 
possible measurement inaccuracy associated with CTA X-wire 
and experimental flow settings [13]. This results an 
overestimation of the mean velocity components than their 
corresponding true values. The deviations may also be partly 
attributed to the inability of RANS approach to accurately 
capture highly swirling flows. 



(a) 



(c) 

Fig. 4 Different turbulence models against experimental data 
[4] for (a) U/U_b for Q_r= 0, (b) U/U_b for Q_r= 1 and (c) 
W/U_b for Q_r= 1. 

3. Results and Discussion 

This section includes presentation of mean and turbulence 
controlling parameters, such as velocity, boundary layer 
thickness, pressure drop, wall shear stress and turbulent shear 
stress for the effect of different number of tangential ports at 
two flow conditions (Q r = 0 and Q r = 1). 

Swirl decay along the length is shown in Fig. 5 for the 3TP 
case. The swirl number is calculated using the Equation (7) and 
Equation (8). It is observed that near the nozzle inlet the swirl 
number S is highly fluctuating and after the converging section, 
the swirl number is nearly constant. Initially, the swirl number 
decreases from the inlet. When the tangential port has 
introduced the value of swirl number increases and 
immediately after the increment the swirl number drops again. 
In the nozzle converging section the magnitude of the swirl 
number rises again at first and then the value decreases and 
finally increases before coming at a constant magnitude. In case 
of S' the values initially increasing from the nozzle inlet come 
to a constant magnitude after the nozzle converging section. It 
is evident that the swirl number S shows a very unpredictable 
nature before the nozzle converging section. 
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Fig. 5 Swirl decay along the length of the nozzle for 3TP case. 

Fig. 6 shows the local swirl numbers along the radial 
location of the nozzle at the different axial positions ( x/D = 
1.275, 3.1, 5.65,9.4, and 14.425) of the nozzle for three 
different numbers of tangential ports, namely, 2TP, 3TP, and 
4TP. It is observed that the profile of local swirl at x/D = 
5.65,9.4, and 14.425 are almost the same in nature for all the 
cases. The local swirl increases gradually from the center of the 
nozzle toward the wall and near the wall, a sudden drop occurs. 
However, the local swirl profiles at x/D = 3.1 are similar for 
2TP and 3TP although for 4TP the local swirl profile is 
parabolic from the nozzle center toward the nozzle wall. The 
value of local swirl near the inlet at x/D = 1.275 is highly 
unpredictable due to the swirl induction. The local swirl profile 
for 2TP and 3TP are nearly similar but for 4TP the profile is 
different. The local swirl near the nozzle wall suddenly 
increases and drops immediately after the increase and finally 
raise to the wall for both 2TP and 3TP. However, the local swirl 
profile at 4TP is like a bell shape; increasing from the nozzle 
center it becomes constant at r/D = 0.2 to r/D = 0.35 and 
then decreases again towards the nozzle wall. As the swirl is 
introduced at this location the local swirl profile becomes 
highly unpredictable. 

The axial velocity vectors at five different axial locations 
(xD = 4.75, 7.50, 9.75, 12.5 and 14.425) are presented in Fig. 
7 for non-swirling (Fig. 7a) and swirling (Fig. 7b-d) 
conditions. Velocity is found to be the almost uniform after the 
converging section for no swirl flow (Q r = 0), but the 
magnitude is the highest near the exit plane, with a reduction 
towards the wall due to boundary layer formation. In contrast, 
for swirl flow (jQ r = 1), velocity magnitude is found to be the 



higher near the inlet and towards the nozzle exit, the velocity 
magnitude decreases. The axial velocity magnitude is 
proportional to the number of tangential ports and the velocity 
magnitude is zero at the wall due to no-slip condition. The axial 
velocity vectors at x/D = 4.75 is symmetric from the 
centerline when swirl is induced (Fig. 7b-d), however, the 
velocity distribution becomes asymmetrical as the flow 
approaches the outlet. The velocity distribution is more 
symmetrical along the radial direction for 2TP than 3TP & 4TP 
with a slight decrease near the center except for exit plane. At 
x/D = 14.425, velocity reduction near the center is highest for 
2TP. 



0 0.1 0.2 0.3 0.4 0.5 

r/D 


Fig. 6 Local swirl charactersitics along the length of the 
nozzle (a) 2TP, (b) 3TP, (c) 4TP. 
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Fig. 7 Axial velocity vectors at different axial locations for non-swirling (a), and swirling flows: (b) 2TP, (c) 3TP and (d) 4TP. 
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Fig. 8 depicts the boundary layer thicknesses along the 
length of the nozzle for both non-swirlign and swirling (2TP, 
3TP and 4TP) conditions. It is observed that the magnitude of 
boundary layer thickness at 2TP is pretty low compare to the 
other. For swirl flow conditions the boundary layer thickness at 
the beginning is small and near the exit the boundary layer 
thickness increases. At the middle of the nozzle the boundary 
layer thickness is nearly constant. However, for no swirl flow 
condition the boundary layer thickness higher and then the 
boundary layer thickness gradually decreases towards nozzle 
exit. 



3 6 9 xW 12 15 


Fig. 8 Boundary layer thickness along the length of the 
nozzle. 

Fig. 9 demonstrates the axial ( U/Ub ) and tangential ( W/Ub ) 
velocity components near the wall i.e. at y/D = 0.01,0.05 
and 0.1 for Q r = 0 and Q r = 1. In the figure, Y represents the 
distance from the wall of the swirl nozzle. It is observed that at 
for swirl condition the near-wall axial velocity profile is 
identical for all the location and it increases from the inlet to 
the nozzle converging section. Then, the velocity suddenly 
decreases at Y = 0.0ID immediately after the converging 
section. The axial velocity is then nearly constant for all the 
axial location of the nozzle. However, the axial velocity 
decreases very little after the nozzle convergence at Y = 0.05 D 
and 0.1D, and no significant change occurs. For swirl flow 
condition the near-wall axial and tangential velocity increases 
after the swirl flow introduction. Then after the nozzle 
converging section the near-wall axial and tangential velocity 
show a very little deviation and along the nozzle axial location 
the velocity profile is nearly uniform at Y = 0.05 D and Y = 
0.1 D. However, the near-wall axial and tangential velocity 
profile at Y = 0.0 ID is fluctuates highly along the length of the 
nozzle, especially the tangential velocity profile shows a totally 
unpredictable nature. The overall near-wall axial and tangential 
velocity distribution at Y = 0.0ID is low in magnitude than at 
Y = 0.05D and Y = 0.1 D. It can be observed that the velocity 
profile at the wall vicinity (Y = 0.0ID) response highly than 
the other position. 

The interplay between the number of tangential ports and 
the swirl intensity for the effect of skin friction coefficient 
along the length of the nozzle is shown in Fig. 10. It is observed 
that for no swirl condition the magnitude of Skin friction 
coefficient is very low and along the nozzle axial location no 
change occurs. However, the Skin friction coefficient at high 
swirl flow at 2TP is the maximum along the nozzle axial 
location although for 3TP and 4TP no significant deviation is 


observed. In every case of swirl flow the Skin friction 
coefficient at the nozzle exit is the maximum. Although the 
Skin friction coefficient along the axial location is not 
significantly changing however, near the exit it suddenly 
increase. The overall distribution of the Skin friction 
coefficient is not similar for all the case. 

Fig. 11 displays the pressure drop at different radial 
locations for swirling flow (3TP). The inlet pressure is taken at 
the reference pressure and pressure drop is calculated based on 
the pressure of the nozzle inlet. It is observed that the pressure 



Fig. 9 Near-wall normalized axial (U/Ub) and tangential 
(W/Ub) velocity distrution along the nozzle length: Q r = 0 (a), 
and Qr— 1 (b and c). 

drop at the center of the nozzle (r/D = 0)is greater than all 
other radial location in almost every axial position of the 
nozzle, however, near the nozzle outlet it suddenly reduces. 
While pressure drop at every radial location increase after the 
introduction of the tangential ports, at (r/D = 0.4) it suddenly 
decrease. Pressure drop near the wall is gradually increasing 
after the converging section of nozzle but pressure drop near 
the nozzle center is nearly constant. It can be concluded that 
Due to the tangential flow introduction the pressure drop 
characteristics near the nozzle center response quickly and near 
the wall vicinity this property change slowly. 
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Fig. 10 Skin friction coefficient along the nozzle length. 



Fig. 11 Pressure drop along the nozzle length for swirling 
flow. 



U1 U4 ^ 0 0.1 0.2 o .. 

r/D r/D 

Fig. 12 Turbulent normal stress profiles at different axial locations. 


Fig. 12 presents the radial distribution of normalized 
turbulent normal stress components at various axial locations. 
It is observed that the turbulent flow characteristic near the wall 
is very high at x/D = 9.4 and x/D = 14.425 at no swirl flow 
( Q r = 0) although at the center of the nozzle the turbulent flow 
characteristic is very low. Moreover, near the inlet of the nozzle 


at x/D = 1.275 and x/D = 4.45 the overall turbulent flow 
characteristic is very small in magnitude. On the contrary at 
swirl flow condition ( Q r = 1) near the inlet of the nozzle at 
x/D = 1.275 the overall turbulent flow characteristic is very 
high. However, the turbulent normal stress near the outlet 
section at x/D = 9.4 and x/D = 14.425 is small in these 
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locations. While at no swirl condition the turbulent normal 
stress characteristic at the center of nozzle is very low at swirl 
flow condition, the u'u' and w'w' component is very high at 
the nozzle center especially at x/D = 1.275 due to the 

<?r = 0 


introduction of swirl flow. The overall magnitude of the 
turbulent normal stress is very small in no swirl flow condition 
but for swirl condition it is very high. 

Qr = 1 





r/D r/D 


Fig. 13 Turbulent shear stress profiles at different axial locations. 


Finally, Fig. 13 presents the radial distribution of 
normalized turbulent shear stress components at various axial 
locations. It is observed that for no swirl flow ( Q r = 0) the 
turbulent shear stress component u'v' and l/w'near the wall is 
very high at x/D = 9.4 and x/D = 14.425 although at the 
center of the nozzle the turbulent share stress component is very 
low (Fig. 13a and Fig. 13c). However, the u'w' component at 
x/D = 4.45,9.4 and 14.425 is reducing near the wall of the 
nozzle. Near the inlet of the nozzle at x/D = 1.275 the overall 
turbulent share stress is uniform and does not deviate much for 
no swirl flow ( Q r = 0). On the other hand at swirl flow 
condition ( Q r = 1) near the inlet of the nozzle at x/D = 1.275 
all the three component of turbulent share stress is fluctuates 
highly and the magnitude of overall turbulent share stress is 


very high. However, the overall turbulent share stress at x/D = 
9.4 and x/D = 14.425 is very uniform along the radial 
location for swirl flow. The overall magnitude of the turbulent 
flow characteristic is very small in for swirl flow compared to 
the swirl condition. 

4. Conclusion 

An incompressible turbulent swirling air jet is investigated 
numerically. In this regard, an axial-plus-tangential swirl flow 
is considered for non-swirling (Q r = 0) and highly swirling 
(Q r = 1) cases. Finite volume based commercial software 
ANSYS Fluent vl7 is used in the simulation to investigate 
mean flow and turbulence characteristics. Governing equations 
are approximated by the RANS equations and turbulence is 
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characterized by the SST k-co model. The study examines the 
effect of number of tangential inlets on mean flow behaviors 
and turbulent characteristics. The magnitude of swirl decay is 
fluctuates highly before the nozzle converging section and after 
the nozzle converging section the swirl decay is nearly 
constant. The boundary layer thickness in swirl flow is found 
to be the smallest for 2TP. The skin friction coefficient along 
the nozzle axial position at no swirl flow is uniform and the 
magnitude is very small. For swirl flow, the skin friction 
coefficient along the nozzle axial position is fluctuating 
especially for 2TP. Moreover, the skin friction coefficient at the 
nozzle exit is the maximum for all the case. Due to the 
tangential flow introduction, the pressure drop near the nozzle 
center response quickly and near the wall vicinity this property 
changes slowly. The turbulent normal and shear stress for no 
swirl flow is nearly uniform, but for swirl flow it fluctuates the 
most near the section where the swirl flow is introduced. 
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ABSTRACT 


This paper evaluates the thermal hydraulic behavior of a pressurized water reactor (PWR) when subjected to the event of Loss of 
Coolant Accident (LOCA) in any channel surrounding the core. The accidental break in a nuclear reactor may occur to circulation pipe 
in the main coolant system in a form of small fracture or equivalent double-ended rupture of largest pipe connected to primary circuit 
line resulting potential threat to other systems, causing pressure difference between internal parts, unwanted core shut down, explosion 
and radioactivity release into environment. In this computational study, LOCA for generation III+ VVER-1200 reactor has been carried 
out for arbitrary break at cold leg section with and without Emergency Core Cooling System (ECCS). PCTRAN, a thermal hydraulic 
model-based software developed using real data and computational approach incorporating reactor physics and control system was 
employed in this study. The software enables to test the consequences related to reactor core operations by monitoring different 
operating variables in the system control bar. Two types of analysis were performed -500% area break at cold leg pipe due to small 
break LOCA caused by malfunction of the system with and without availability of ECCS. Thermal hydraulic parameters like, coolant 
dynamics, heat transfer, reactor pressure, critical heat flux, temperature distribution in different sections of reactor core have also been 
investigated in the simulation. The flow in the reactor cooling system, steam generators steam with feed-water flow, coolant steam flow 
through leak level of water in different section, power distribution in core and turbine were plotted to analyze their behavior during the 
operations. The simulation showed that, LOCA with unavailability of Emergency Core Cooling System (ECCS) resulted in core 
meltdown and release of radioactivity after a specific time. 


Keywords: LOCA; VVER-1200; PCTRAN; ECCS; Thermal-hydraulic behavior. 



This work is licensed under a Creative Commons Attribution-NonCommercial 4.0 International 


1. Introduction 

VVER-1200 is the utmost type of VVER pressurizer water 
reactor (PWR) meaning Vodo Vodyanoi Energetichesky 
Reaktor or water-water energetic reactor with improved safety 
features design and developed by OKB Gidropress, Russia. It is 
a thermal neutron pressurizer reactor used water not only as 
coolant but also as moderator. It consists of two circuit namely 
the primary circuit, where the primary coolant receives heat from 
reactor core and the secondary circuit where the heat supplied to 
produce electricity [1] with four cooling loops, main circulation 
pump, accumulator tanks of the ECCS, relief and safety 
emergency valves. It’s said to be a generation III+ reactor with 
improved safety features like sacrificial concrete core catcher 
layer, passive heat removal system, etc. Inside the core, the 
primary coolant deposit at 155 bar so that it can avoid 
vaporization at high temperature [2]. Any leakage in the coolant 
circuit can lose large amount of coolant from Reactor Cooling 
System (RCS) result in Design Basis Accident or Beyond Design 
Basis Accident according to the flow from RCS. 

Design-basis accident (DBA) is a postulated accident that a 
nuclear facility must be designed and built to withstand without 
loss to the systems, structure and components necessary to ensure 
public health and safety. DBA concept is very feasible for 
abnormal operations in a nuclear power plant to vindicate 
performance requirements for reactor system, structure and 
components after anticipated operational occurrence. Beyond 


Design Basis Accident (BDBA) which have extremely low 
probability of happening is a sequence of accident which may or 
may not lead to core degradations. DBA postulate for each type 
of reactor covering different type of failure combinations like 
reactivity-initiated accident, loss of flow transients, Loss of 
coolant accidents, fuel handling accidents and so on. The DBA 
is prescribed to different limiting condition for cladding 
temperature, local fuel cladding oxidation, ZrC >2 mass cladding, 
whole body dose and others [3]. 

LOCA is one of most important and common DBAs for 
each reactor design depending on size of break. LOCA assumed 
that one of inlet or outlet pipe from circulating pump to reactor 
vessel completely broken or partially broken and free discharge 
of primary coolant from both broken ends occur. LOCA can be 
happening for different reasons like, material defects of different 
systems, material fatigue in pipe/shell, external impacts like 
heavy loads or missiles, device failure during operation due to 
different reasons, inadvertent opening on primary system 
boundary and so on. The nuclear reactor core has different 
emergency system to mitigate different type of accident depend 
on the incident or accident type. 

There are various safety features in VVER 1200 reactor 
plants. VVER 1200 reactors are capable of sustaining initiating 
events as well as their consequence under accident conditions 
during a long period within the boundaries of design safety 
criteria and assured by following structural and design features 
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like (a) increased coolant volume above the core (b) increased 
coolant volume in the primary circuit in respect to fuel mass and 
thermal power of the core (c) increased capacity of pressurizer 
(d) reliable natural circulation (e) considerable water inventory 
in horizontal SGs on the secondary side [2]. 

The primary coolant system of reactor core is the system 
which is heated at reactor core flows into steam generator. The 
primary coolant flow through heat exchanger tube boiling the 
coolant in secondary cooling system which flow outside of heat 
exchanger tube that drive the turbine generator and creates 
mechanical energy from heat energy. After that, the primary 
coolant comes from steam generator flow back to reactor core 
using compressed pump. 

The leakage size within reactor coolant pressure boundary 
divided into two group namely Large Break Loss of Coolant 
Accident (LBLOCA) and Small Break Loss of Coolant Accident 
(SBLOCA). Reduction of water temperature during Emergency 
core cooling system injection drop in about 50°C-70°C the break 
can be categorized as small breaks and for large break reduction 
of temperature about 120°C-150°C [4]. The primary system 
remains at high pressure after the occurrence of LOCA and the 
core remain covered. Coolant flows inside the primary system as 
soon as the pumps are tripped due to automatically or manually 
gravity control phase separation occurrence. The next 
subsequent events depend on the overall behavior of the primary 
and secondary system. Automatic and operator-initiated 
mitigation measurements influenced the overall behavior. 
Reactor design, break size, location of the break, core bypass 
system, operator interactions have influenced the scenario. Each 
alternative action has much impact to make a huge difference on 
the LOCA. The LBLOCA classical DBA for reactor core where 
one of the inlet pipes from the circulating pump is completely 
non-functional or free discharge of primary coolant from both 
broken ends which referred as Double Ended Guillotine break. 
LOCA have those following phase (a)BLOWDOWN PHASE, 
0-20 seconds (b) BYPASS PHASE, 20-30 seconds (c) REFILL 
PHASE, 30-40 seconds (d)REFLOODING PHASE, 40-250 
seconds (e) LONG TERM COOLING >250 seconds. 

The VVER-1200 has following active and passive system to 
mitigate different accident: (i) The Primary Circuit Emergency 
with Planned Cooldown and Spent Fuel Cooling (ii) Low 
Pressure Emergency Injection System to supply Boric Acid 
Solution (iii) Emergency Core Cooling System (iv) Passive Core 
Flooding System (v) Emergency Boron Injection System (vi) 
Emergency Gas Removal System (vii) Primary Overpressure 
Protection (viii) Secondary Overpressure Protection (ix) Passive 
Heat Removal System (x) Steam Generator Emergency 
Cooldown System (xi) Main Steam-line Isolation System (xii) 
Double Envelope Containment and Core Catcher [1]. 

Personal Computer Transient Analyzer (PCTRAN) is a 
nuclear reactor transient and accident simulation tool that 
operates on a personal computer for learning different nuclear 
reactor behavior under different condition. Concepts of neutron 
multiplication, criticality, thermal hydraulic safety parameters, 
transient analysis, accidental dose estimation and Xenon 
poisoning and feedback in live simulation statistical analyses can 
be carried out. The PCTRAN-VVER1200 simulator version 
1.2.0 available from Micro-Simulation Technology for VVER- 
1200 reactor for thermal output of 3200 MWt with two loops of 


reactor coolant system represented as loop A and loop B use for 
the numerical simulation purpose. It has shown the major 
systems and components of the reactor like RCS, the Emergency 
Core Cooling System (ECCS), the Pressurizer with its 
components, the Low-Pressure Injection System (LPIS), 
Feed/bleed water system, Steam Generator, Accumulator, 
different safety valves and different components of the NPP. The 
components like pump, valves indicating red colored are 
operating and open system and white indicating colored indicate 
idle and closed in the mimic interactive. Those components can 
be switch on or off just by clicking on them in the interactive 
display. It has two interface, the main interface for different 
component and Dose mimic interface for radioactivity display. 

Transient behavior simulation has been successfully 
benchmarked by using PCTRAN approved by International 
Atomic Energy Agency (IAEA). PCTRAN based on end-users 
input data with a reactor point kinetics model for analysis nuclear 
safety related concepts and phenomena that could be easily 
explained with the simulators along with understanding of 
valuable technological differences in various designs of 
PCTRAN products. A high-resolution color mimic of the 
Nuclear Steam Supply System (NSSS) and containment displays 
are the status of important parameters and allows simulation of 
operator actions by interactive control. 

Though the main purpose of this simulator is operator 
training and a dynamic test to validate the control logics in 
reactor regulating system (RRS), this simulator is also used for 
research purposes. Like, Mollah review about PCTRAN 
simulator for safety and transient analysis of PWR [5]. Ibrahim 
et al. showed Simulation of Safety and Transient Analysis of a 
Pressurized Water Reactor using the PCTRAN [6]. Yi-Hsiang 
Cheng et al. integrated development of accident dose 
consequences simulation software for nuclear emergency 
response applications using PCTRAN software [7]. Also further 
they analyze the Improvement of accident dose consequences for 
nuclear emergency response applications using PCTRAN 
simulator [8] . Khan and Islam published an article on a PCTRAN 
based investigation on the effect of inadvertent control rod 
withdrawal on the thermal-hydraulic parameters of a VVER- 
1200 nuclear power reactor [9]. Fyza et al. analyzed the thermal- 
hydraulic parameters of VVER-1200 due to loss of coolant 
accident concurrent with loss of offsite power with PCTRAN and 
conclude with PSAR (Preliminary Safety Assessment Report) 
regarding LOCA [10]. Tube rupture in steam generator and 
transient analysis of VVER-1200 using PCTRAN simulation is 
conducted by Saha et al. [11]. 

Only a limited research works on VVER-1200 new 
generation III+ reactor related topic has been conducted. In 
this study, the thermal hydraulics behavior of VVER-1200 
reactor due to LOCA with 500% area with and without ECCS 
are simulated by using PCTRAN software. Thermal hydraulic 
parameters like, coolant dynamics, heat transfer, reactor 
pressure, critical heat flux, temperature distribution in different 
sections of reactor core have also been investigated in the 
simulation. The flow in the reactor cooling system, steam 
generators steam with feed-water flow, coolant steam flow 
through leak level of water in different section, power 
distribution in core and turbine were plotted to analyze their 
behavior during the operations. 


54 




P.D. Nath et al. /JEA Vol. 01(02) 2020, pp 53-60 


2. Methodology 

The numerical analysis was carried out by PCTRAN 
VVER1200 v-1.2.0. By this software package one can simulate 
a variety of accident and transient condition for VVER-1200 
reactor nuclear power plants. User can select initial conditions 
corresponding to different parameter and initiate malfunctions 
to analyze as desired. Following analysis had been conducted 
by this software- normal operating control, stream line break or 
LOCA, recirculation pump trip, turbine trip, station blackout, 
inadvertent rod withdrawal or insertion, Boron dilution 
transient, steam generator tube rupture, feedwater transient, 
anticipated transient without scram and any combination of 
above. The user can manually trip the reactor pump, open or 
close valve, override ECCS, change operational set points. The 
plant mimic system has the ability to freeze, back track, snap a 
new initial condition, change the simulation speed, trend plot 
by selected variable for conducting system analysis. 

The simulator based on nuclear reactor theory and solution 
technique by Six- group point kinetics model which is a point 
kinetics equation with six delayed neutron groups and 
reactivity control from external sources and feedback. The 
model equations are described in the literature by A.S. Mollah 
[5]. The initial steady state conditions of the plant were Reactor 
Power at 100% with RC Pressure 155 bar and Core Average 
Temperature 306.9°C with BOC (Beginning of Cycle) of Steam 
Generator (SG) pressure 70 Bar. After setting up the initial and 
boundary conditions, the necessary operating parameters were 
measured for the analysis of reactor performance. Default 
values of VVER-1200 reactor were given which can be 
changed to initial value as required. Different type of 
malfunctions with initial value, fraction percentage, delay time 
will be set initially to run simulation. 

After completion of simulation, transient plot of different 
variables like flow rate, coolant activity, concentration of 
materials, dose rate, fraction, energy, level of water, mass of 
different element, power of different system, pressure, 
radioactivity, enthalpy temperature, void fraction could be 
created with respect to time and other variables. The sequence 
of simulation of the whole process is shown in Fig. 1. 

Step 1. Set and active the simulation for run 

Step 2. Default initial condition or set new 
initial condition as required 

Step 3. Describe malfunction with delay 
time and failure fraction 

Step 4. Run the simulation in control panel 

Step 5. View transient plot and data for 
different variables 

Step 6. Save transient report and dose data 
for analysis 


Fig. 1 Steps/order of simulation in PCTRAN. 


3. Computational Results and Discussion 

In this computational analysis two cases have been 
considered for investigation. First case accounted for a design 
basis loss of coolant accident where the ECCS was in operation 
and the second case was for a beyond design basis loss of 
coolant accident without any type of ECCS and passive cooling 
system. After conducting the computational analysis, the plot 
relevant to different parameters were discussed. 

3.1 Analysis of LOCA at cold leg with all ECCS 

In this simulation 500 cm 2 (25.23 cm diameter) break in 
cold leg simulation have been conducted. LOCA with break at 
cold leg has serious impact because safety injection system 
makeup water lost through this break. The simulation carried 
out using manually defined malfunction in the code control at 
cold leg with 500 cm 2 failure fractions after 5 seconds of the 
simulator to run at steady state condition as shown in Fig. 2. To 
minimize the coolant loss, the RCPs manually tripped. There 
was a rupture icon in the cold leg with primary coolant loss rate 
appears in the simulator. Safety injection systems started 
automatically after reactor trip signal generated on because of 
the low pressure. 



Fig. 2 Graphical user interface of PCTRAN: Initiation of 500 
cm 2 cold leg LOCA. 



Fig. 3 Graphical user interface of PCTRAN: Reactor trip and 
actuation of HPIS, RB spray pump. 
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Fig. 4 (a) Pressure of RCS, RB, SGs (bar) for cold leg break (b) Mass flow rate of RCS leak, reactor coolant loops, ECCS flow 
water (kg/s) (c) Temperature of average fuel, RCS, peak clad, peak fuel, RB (°C) for cold leg break (d) Core thermal power (%), 
Nuclear flux power (%), Turbine Load power (%) for cold leg break (e) Water level of core, pressurizer, RB sump, SGs (m) for 

cold leg break. 
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Fig. 3 shows that, ECCS actuated High Pressure Injection 
System (HPIS) at 15s when the RCS pressure was lower than 
the nitrogen pressure accumulator tank. Although HPIS along 
was not enough to compensate primary coolant loss, 
accumulator and Low-Pressure Injection System (LPIS) with 
recovery refueling water storage tank, reactor building spray 
help to provide boronated water. LPIS provide enough water to 
maintain the core submerged and responsible for long term 
cooling of residual heat in reactor. Compared to large break 
LOCA, the depressurization occurs relatively slowly in 
SBLOCA. Upon actuation of accumulator, ECCS water flow 
started to exceed the coolant discharge rate and the reactor 
vessel water level slowly rose again and core level water 
restored as illustrated in Fig. 4 (e). The RCS pressure fall down 
from initial steady 160 bar to about 3 bar after 200 second from 
the onset of accident. The Steam Generator A (SG-A) and SG- 
B pressure initially increased after accident from 74 bar to 81 
bar in 30 second for sudden loss of coolant but after 65 second 
both SGs reached to equivalent pressure 72 bar as shown in Fig. 
4 (a). The process of depressurization was very rapid slowing 
down the reactor coolant leak rate. After being emptied the 
pressurizer during initial occurrence, RCS leak rapidly 
decreased as the pressure altitude change fast. RCS leak 
initially started after 5 seconds of the incident with 4150 kg/s 
mass flow rate which slowly decreased with respect to time 
after 40 seconds of the accident. The leak loss of mass flow 
become equal with reactor coolant loop B flow at 105 seconds 
which was about 390 kg/s. After 220 seconds the mass flow 
rate become very slow but it needed more time to stop 
completely. ECCS started after 85 seconds with gradual mass 
flow rate highest of 1390 kg/s. After 245 seconds, the ECCS 
slow down. The reactor coolant loop A has 4300 kg/s mass flow 
rate, coolant loop B has highest 13000 kg/s mass flow rate 
which decrease very rapidly and stable in 150 seconds to 225 
seconds at 390 kg/s in Fig. 4 (b). For SBLOCA initiation 
reactor coolant release very high flow rate through break, 
which loss rate decrease due to HPIS actuation and 
depressurization of primary system. 

In the meantime, emergency accumulator water system 
active to provide enough makeup water in RCS to maintain the 
core submerged in water so that during the transient event, fuel 
and cladding temperature never rise to critical value shown in 
Fig. 4 (c). Highest peak clad temperature and peak fuel 
temperature can’t rise to critical value to melt the core as the 
control rods are fall down after 7 s and huge makeup water from 
ECCS, HPIS active due to this transient event. Fuel peak 
temperature, clad peak temperature, average fuel temperature 
falls from respectively 800°C, 415°C, 320°C to 150°C, 140°C, 
110°C after 300s. The sump water level in reactor building 
increase from 40°C to 100°C due to recirculation through 
RWST with LPIS. The accident is assumed to happened when 
the reactor operating at full 100% power. Due to the occurrence 
of LOCA in cold break, the simulation shows reactor trip after 
40s which indicate core thermal power and nuclear flux fall to 
5% and decay heat power continuously extend for next few 
cycles to cooled down to zero in Fig. 4 (d). The turbine load 
initially dropped after LOCA which increase for following 20s 
and again decrease to zero as turbine trip after the malfunction. 


3.2 Analysis of LOCA at cold leg without ECCS or any passive 
heat removal system 

In VVER-1200, the emergency core cooling system is 
consisting with not only HPIS, but also it has LPIS, passive heat 
removal system, emergency accumulator water system, feed 
water system, reactor building spray pumps and refueling water 
storage tank with it to prevent all kind of accident including 
DBA and BDBA to tackle the situation in emergencies. 
Whenever ECCS face any kind of accident signal, it 
automatically got actuated to keep the core under water and 
cooled it down in order to eliminate any serious impact on the 
total system and environment with working personnel. When 
the ECCS can’t operate as needed, then this scenario can turn 
into severe accident with core meltdown and fuel degradation. 

In this case, all the ECCS equipment were disabled and 
reactor coolant pump was turned off manually, with a 500 cm 2 
LOCA in cold leg. Though, it is a hypothetical cause with 
unavailability of all those ECCS equipment but following the 
Fukushima Daiichi accident caused by earthquake and tsunami 
where all the emergency system failed and level/class 7 disaster 
happened, possible consequences of this type of accident needs 
to be analyzed. 

For LOCA accident without enabling the ECCS, the 
reactor trips instantaneously and after 9 seconds all the control 
rods trip down to reactor core with pressurizer 
emptied within 20 seconds due to very rapid depressurization 
as shown in Fig. 5 (a). The fuel temperature continuously 
increases due to absent of all ECCS and passive heat removal 
system. With increased cladding temperature, the interaction 
between the cladding and steam was accelerated. Hydrogen 
concentration inside reactor increased with time and shortly 
after the core uncovery, as depicted in Fig. 5 (b) about 980 
seconds after the accident, the hydrogen concentration in the 
reactor building exceeded 5%. Fuel temperature increased as 
time goes by resulting the reactor core starts to melt. In Dose 
Mimic, the core melting process can be observed in Fig. 5 (c). 
Molten fuel collapse into the bottom of the vessel as there is 
nothing that can cooled the core materials. The vessel lower 
head may then heat up to the melting point, too. Fig. 5 (d) 
illustrates that, the molten debris with molten core concrete 
interaction with vessel penetration may drop into the 
containment cavity floor takes around 3000 seconds for 500cm 2 
LOCA at cold leg without ECCS. During the fuel damage 
process, first the fission gas in the cladding may leak out. Later, 
if the fuel and cladding continue their degradation, 
radionuclides will also be released. In addition to iodine and 
noble gases, there are alkali metals, tellurium, barium, cerium, 
lanthanides, etc. The elevated concentration of these 
radionuclides would find their ways through the vessel break, 
relief valves, and containment leakage into the environment. 
The molten metal interacts with concrete and forms a slump. At 
lower temperatures, degassing of concrete occurs and both 
steam and carbon dioxide can be released. At higher 
temperatures concrete can also be molten and mixed with 
metals. In the event that the molten core heats up the vessel 
bottom and melts through it, the debris, called corium, falls into 
the reactor cavity. 
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(c) (d) 

Fig. 5 Graphical user interface of PCTRAN: (a) Reactor trip for LOCA without ECCS (b) Complete core uncovery (c) Beginning 
of core melt (d) Molten core-concrete interaction (MCCI) with vessel penetration. 


For LOCA at cold break, the reactor coolant system leak 
and loops fluid decrease rapidly. At the time of emptying of 
pressurizer in first few seconds, coolant leak rate rapidly 
decreased as the pressure elevation of break firstly change. 
From Fig. 6 (a), level of reactor building sump slowly grew up 
to a certain point until the recirculation process terminated and 
there left no water to carry the poison. Level of core water 
started to decrease very slowly at about 73 seconds. There 
created a tube rupture loop at cold leg after 350 seconds. In Fig. 
6 (b), different temperatures were plotted as respect to time for 
initial 300 seconds. Due to increasing the sump in reactor 
building, temperature also increased as followed by the sump. 
Initially after the accident, the control rod fell down to the core 
but without the activity of any coolant flow inside the core, the 
power slowly rose by chain reaction. Before that, the core water 
level and control rods minimized the temperature about 140 
seconds until all the remaining water vaporized totally. After 
the core fuel emersion, temperature rapidly increased due to 
unavailability of coolant. Increase in Zircaloy cladding 
temperature with fuel temperature accelerated the reaction 
between Zirconium and water. Hydrogen gas accumulation 
inside reactor building occurred quickly at the top of building 
which can get exploded when it increased to 5%. The 
Zirconium cladding distortion increased very high; fraction rate 
about 42 in 1500 seconds whereas Hydrogen buildup was 8% 


at that time. DNBR or departure from nuclear boiling is the 
ratio of critical heat flux at specific location to operating local 
heat flux at that location. In all NPP, DNBR must be greater 
than one but it can’t allow to increase more than 95. From the 
Fig. 6 (c), it is observed that the DNBR increased very slowly 
as local heat flux can’t able to cross limit until all the coolant 
system loss and its critical value was about 80. When the local 
flux was higher than critical heat flux, the core materials melted 
down for increasing heat flux. Fuel Doppler reactivity is the 
fractional change in reactivity caused by fuel temperature. Until 
it’s less than one, there was no harm of reactivity but after 
explosion of the reactor building beyond 980 seconds, reactive 
substance released to atmosphere. As it increased after the loss 
of coolant accident, it progressed very much and spread 
radioactive materials as shown in Fig. 6 (d). Fig. 6 (e) shows 
the comparison between LOCA for different temperature in 
fuel (TF), average fuel (TAVG), peak fuel (TFPK), peak clad 
(TPCT) with and without ECCS for first 300s. Initially after the 
accident, all ECCS activated when the control rod submerged 
with fuel inside the core, hence temperature fell down for both 
conditions. But unavailability of ECCS in second case, 
temperature began to rise when all the water inside core 
vaporized and fuel with control rod became unstable due to 
active heat generation. 
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Fig. 6 (a) Level of core, pressurizer, SGs, RB sump water for LOCA at cold leg without ECCS (b) Temperature of average fuel, 
RCS, peak clad, peak fuel, RB (°C) for LOCA cold leg break (c) DNBR for LOCA without ECCS (d) Fuel Doppler (%), rod (%), 
modifier reactivity (%) for LOCA without ECCS (e) Temperature of different component with and without ECCS for LOCA at 

cold leg. 


59 

















































































P.D. Nath et al. /JEA Vol. 01(02) 2020, pp 53-60 


4. Conclusion 

In this study, numerical study has been conducted to 
investigate the different thermal hydraulics characteristics of 
nuclear power plant due to loss of coolant accident with and 
without ECCS. Loss of coolant accident may be a type of 
postulate event but the probability of happening it in any kind 
of nuclear reactor plant have enough chance due to different 
reasons. LOCA accident with both design basis accident 
analysis and beyond design basis accident analysis are 
employed in this paper. From the results following conclusion 
could be drawn: 

• Design Basis Accident and Beyond Design Basis Accident 
both must be mitigated before they are harmful for 
environment and mankind. A simple Design Basis 
Accident can turn into serious disaster if it won’t handle 
properly. Risk of radioactive material can be spared if 
proper maintenance isn’t taken care off. Serious accident 
is not only dangerous for human, but also for every life on 
that area. 

• Different type of thermal hydraulics behaviors like 
pressure, flow rate, water level in different sections, clad 
fuel temperature, feed water cooling system, power and 
many things can be investigated through LOCA analysis. 
From numerical simulation for LOCA without ECCS 
transient behaviors and core degradation dose mimic it is 
found that it takes 3066 seconds to complete core 
meltdown for 500% area or (25.23 cm diameter) failure 
fracture. Different arbitrary failure fraction takes different 
time to conclude in MCCI but reactivity release must be 
controlled. Design basis LOCA show good agreements 
with Primary Safety Analysis Safety report (PSAR) but in 
case of BDB LOCA occurring without ECCS, the core 
melting process is not able to mitigate the accident. 
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ABSTRACT 


Cancer is one of the most leading causes of death now worldwide. Most of the cancer therapy aims to raise the temperature of the 
cancerous tissue above a therapeutic value and thermally kill or destroy it. Minimizing the damage of the healthy cells surrounding the 
infected cells is one of the major concerns of these therapies. Precise acknowledgment of the temperature profile of living tissue during 
therapy is of utmost necessity for this purpose. Towards that direction, this paper presents an unsteady finite element model of the 
bioheat equation to analyze the temperature distribution during the thermal therapy. A C language based system has been developed to 
solve the unsteady part of the problem employing Crank-Nicolson method and to solve the linear problem employing the Gauss 
elimination technique. Using this system, we investigate thermal behaviors in human tissues subjected to constant, sinusoidal spatial 
and surface, point, and stochastic heating. It was found that surface heating is beneficial for treating skin surface cells, while laser 
heating for the cells that lie below the skin surface. Moreover, for deep cell, the point heating style can bring the most desirable outcome. 
Results describe in this paper could be useful for researchers and doctors to optimize the treatment procedure, even protocols. 


Keywords: FEM; Bioheat Transfer; Human Tissue; Pennes Equation; Spatial Heating. 
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1. Introduction 

Prior acknowledgment of the temperature distribution of 
thermal therapy could help optimize treatment procedures and 
take necessary precautions for probable danger during therapy. 
Undoubtedly experimental study is most welcome in this field. 
However, as human subject involves, the experimental approach 
often became very difficult to perform and result in a hazardous 
situation. Moreover, it is time-consuming and needs costly 
equipment setup. Hence mathematical modeling is preferred in 
that area instead of clinical trials. But the complexity of modeling 
biological systems and treatment procedures makes it so tough in 
practice even impossible. That’s why computational modeling of 
the biological bodies has received considerable attention in the 
research community in the past decade. Rapid advancement in 
computational technology, which enables better accuracy with 
less computational cost, added a new era to this progress. 
Moreover, due to its simplicity in use and low cost, it is widely 
used in simulation of biomedical problems. With the help of 
computer technology and mathematical model, it is possible to 
calculate and visualize the stationary and transient temperature 
inside biological bodies during thermal therapy. 

Heat transfer in living tissue has become a very interesting 
topic for scientists and engineers because of its broad application 
in bioengineering and designing of medical protocols. Many 
therapeutic applications need the proper thermal description of 
the human body. Not only that, but many medical operations also 
rely on engineering methods to determine the safety and risk 
level involved in many surgeries. At present, mathematical 
modeling of Bioheat transfer is widely used in treating tumors, 
cryosurgery, laser eye surgery, and many other applications. The 


success of hyperthermia treatment much depends on the proper 
knowledge of heat transfer in blood perfused tissue [l]-[2]. 

Over the years, several mathematical models have been 
developed to describe the heat transfer within living biological 
tissues. The most widely used bio-heat model was introduced by 
Pennes in 1948 [3]. Pennes Bio-heat-equation has been widely 
used to approximate the overall temperature distribution in 
tissue. The Pennes model was initially used to predict the 
temperature distribution in the human forearm. Due to its 
simplicity (uniform thermal conductivity, blood perfusion, and 
metabolic heat generation), it was implemented in various 
biological research work such as for cancer hyperthermia, 
cryosurgery etc. Therefore, to obtain a flexible solution that can 
solve similar problems is very desirable. Reports on the 
analytical and numerical solutions of the bio-heat transfer 
problem are found in the literature. In some analytical cases, 
sinusoidal heat flux [4] and sometimes cooling of the skin [5] 
were considered as boundary conditions. In some numerical 
analysis, sinusoidal surface heat flux was used as boundary 
condition [6]. Researches related to the bio thermos- mechanical 
were reviewed in [7]. Monte Carlo method was used to solve the 
multidimensional problem in [8]. The result of boundary element 
method and finite element method for the numerical solution of 
the steady state bio-heat transfer model of the human eye were 
compared in [9]-[10]. The finite element method was used for 
the thermal-magneto static analysis in biological tissues in [11]. 
Some commercial software includes bio-heat transfer functions 
with limited boundary conditions. It is quite difficult to have the 
temperature profile of a particular point or a line within the 
biological tissue with different time intervals using that software. 
The development of a free finite element code of bio-heat 
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equations can meet the purpose. The objective of this research is 
to develop a one dimensional finite element code for the solution 
of both steady and transient bio-heat equation. The popular 
Crank-Nicolson method was used in time discretization of the 
transient analysis. The developed finite element code was used 
to simulate the thermal response of tissue during cancer 
hyperthermia, laser surgery, tissue heating with a hot disk, and 
point heating. Moreover, time-dependent spatial and surface 
heating were incorporated. The effect of surface heating, step 
heating, sinusoidal heating, and point heating was thoroughly 
investigated. The temperature profile for all cases is found with 
valuable information for the physicians and researchers. 

2. Bio Heat Transfer 

For the study of bio-heat transfer in human tissue, the most 
useful one is Pennes equation which can be expressed as 


is equivalent to both the governing differential equation as well 
as a certain type of boundary condition. The simplest form of 
the Eq. (1) is 



BT + q = pc 


dT 

at 


(7) 


Where B = cobPb Cb and q = BT a +Q m +Qr. The weak form 
of the differential equation (applying the Weighted Residual 
method) is derived as 


C f, + k 1 77 - BWT - H + + 

(WQ) Xb = 0 (8) 


Where W is the weighted function, and Q is the secondary 
variable. A linear element is considered is this model whose 
temperature function is expressed as 


+ “bPbCbCTa - T) + Q m + Q r = pc^ 


( 1 ) Th(x) — 2f=i<pf(x)1j e 


(9) 


For steady state case the Eq. 1(a) is reduced as 


Using the linear approximation of Eq. (9) finally a linear system 
was derived of the following form 


+ to b p b c b (T a - T) + Q m +Q r - 0 


(2) [C]{T} + [K]{T} = {q} + {Q} 


( 10 ) 


Where p, c, k are respectively the density, the specific heat, 
and the thermal conductivity of the tissue; pb, c b denote density, 
and specific heat of blood, respectively. The cob is the blood 
perfusion, T a the known arterial temperature, and T (x, t) is the 
unknown tissue temperature. Where Q m is the metabolic heat, 
and Q r (x,t) is the heat source due to spatial heating. 

Let the one dimensional problem of length L, where the 
skin surface is defined at x = 0 and the body core at x = L. The 
constant body core temperature is defined as T c , ho is the 
ambient heat convection coefficient between the skin surface 
and the surrounding air and To is the ambient temperature. At 
the skin surface (x = 0), the thermoregulation between the skin 
and ambient air is governed by the thermal convection between 
air and skin. While as the tissue temperature remains constant 
within a narrow limit, i.e., 2-3 cm, the boundary conditions at 
the body core is considered as temperature boundary 
conditions. Thus the boundary conditions for this particular 1- 
D problem can be written as: 

= _h ° [T ° - T( »] at x = 0 o) 

T = T c at x = L (4) 

In some cases, force convection cooling is applied at the 
skin surface to remove excessive heat from the skin surface, in 
such case the boundary condition at the skin surface is defined 
as 

-k^ = —h f [T f - T(x)] at x = 0. (5) 

In some cases, the boundary condition is time dependent. So 
time dependent boundary conditions can be expressed as 

-k^ = Q(t) at x = 0. (6) 

Here Q(t) is the time dependent heat flux. 

2.1 Finite element discretization 

The first step of the finite element discretization is to 
develop a weak form that is a weighted integral statement and 


Where C is the capacitance matrix, K is heat conductive 
matrix and T is unknown temperature and others are known 
vectors. 

2.2 Time Discretization Scheme 

A simple time integration scheme for the Eq. (10) is 
derived by assuming that C and K are constant. In such case, 
the matrix differential equation can be discretized with 
response to time as 

r rn+l_rpn 

C AT + aKT n+1 + (1 - a)KT n = Q + q (11) 

Where T n+1 and T n are the vectors of unknown nodal 
values at times nAT and (n -I- 1)AT, respectively and a is the 
weighting factor. The a must be chosen in the interval between 
0 and 1. When the value of a is considered 0.5, the process is 
called the popular Crank-Nicolson method. The discretized Eq. 
(11) can be written as: 

+ cxk) T n+1 = [c^-(l-a)K]T n + Q + q (12) 

The Eq. (12) was solved using an iterative procedure. The 
initial temperature is known and then the temperature of the 
next step is calculated from the solution of Eq. (12) through the 
Gauss elimination technique. 

2.3 Boundary Conditions and Input Parameters 

Throughout the study at X=L, the temperature boundary 
condition is used. 

However, depending upon the types of heating boundary 
condition (3), (5) and (6) is used at X=0. 

In section 3.2, 3.3, 3.8, and 3.9 boundary condition (3) is 
used at the skin surface. Heat flux boundary condition Eq. (6) 
is used in sections 3.4 and 3.5 at X=0. Where the force 
convection boundary condition Eq. (5) is used in sections 3.6 
and 3.7. The input parameters used in this study is summarized 
in Table 1 [12]. 
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Table 1 Input Parameters. 


Parameters 

Value 

Thermal conductivity (k) 

0.5 w/m 2 

Convection Coefficient (h 0 ) 

10 w/m 2 

Forced convection coefficient (hf) 

100 w/m 2 

Environmental Temperature (To) 

25 °C 

Temperature of the Artery (T a ) 

37 °C 

Body core temperature (T c ) 

37 °C 

Metabolic heat generation (Q m ) 

33800 w/m 2 

Density of blood (pb) 

1000 kg/m 3 

Density of tissue (p) 

1000 kg/m 3 

Specific heat of blood (cb) 

4200 J/kg.°C 

Specific heat of tissue (c) 

4200 J/kg.°C 

Blood perfusion (cob) 

0.0005 ml/s/ml 


3. Results and Discussion 


For the simple thermal analysis, from the skin surface to 
tissue body is enough to consider. So to avoid the 
computational complexity, a ID tissue of length (L) of 30 mm 
is considered as a computational model. 

3.1 Code Verification 

Fig. 1 shows a comparison between numerical result and 
the analytical result obtained from [12] for the steady state case 
considering the Eq. (3-6). The comparison shows a better 
agreement. The boundary conditions are defined as describe in 
Eq. (3) and Eq. (4). 



Fig. 1 Comparison with analytical solution. 

This figure shows that initially, tissue temperature 
increases along with distance due to metabolism, but after 
attaining the highest value, it decreases towards the body core 
as temperature boundary condition is applied to the body core. 
Here the maximum temperature is 45°C, which located about 
11 mm below the skin surface. 

3.2 Spatial Heating 

Laser and microwave therapy are some of the most widely 
used non-invasive techniques to destroy malignant cells. In this 
section, we aim to know the temperature distribution of human 
tissue during heating by laser or microwave. 


In case of heating by laser and microwave, the heat 
absorption rate can simply be approximated by Beers Law, 
which can be expressed as Q r = r]P 0 (t)e~ T]X in which heat 
flux decays exponentially with respect to distance from the skin 
surface [13]-[15]. Here P a (t) is the time-dependent heating 
power on the skin surface, and rj is the scattering coefficient. 
Since Po(t) and r| vary from one apparatus to another, so it is 
important to know the influence of these parameters on tissue 
temperature. P 0 (t) can be either constant and time-dependent. 
In our study, we have considered both cases. 



(a) 



(b) 

Fig. 2 Temperature distribution at different times; (r|= 200 m' 1 
, (a)P o (t)=250 W/m 2 , (b) P o (t)=250+200cos(0.02t) W/m 2 ). 

Fig. 2 depicts the transient temperature at different times 
when tissues subject to two different spatial heating. Fig. 2 (a) 
shows the case of constant spatial heating and while Fig. 2 (b) 
is for sinusoidal spatial heating. In both cases, at the early stage, 
the tissue temperature increases along with the distance from 
the skin surface due to external heating, but later it decreases 
towards the body core. Moreover, there is an inter-cross 
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between temperature curves at different times in Fig. 2 (b), 
which indicates the oscillation of the temperature inside the 
tissue due to sinusoidal spatial heating. These figures also 
reveal that the temperature within the tissue increases along 
with time and finally reached the steady state condition. Here 
the maximum temperature is about 47°C, and that lies 7 mm 
below the skin surface. 

3.3 Effect of Scattering Co-efficient 

Fig. 3 shows the effect of the scattering coefficient, where 
Fig. 3 (a) shows the result for constant heating, and Fig. 3 (b) 
depicts the result for sinusoidal heating. In both cases, the 
larger coefficient results in a higher temperature. Moreover, in 
the case of sinusoidal spatial heating, the larger coefficient 
returns higher amplitude. Fig. 3 (b) indicates the sinusoidal 
effect. Fig. 3 (a) shows that after about 3000 seconds 
(approximately), tissue temperature begins to stabilize. 



36 0 2000 4000 6000 800010000 

Time(s) 

(a) 



(b) 

Fig. 3 Effect of scattering coefficient on temperature response 
at skin surface (f (t) =0); (a) P 0 (t) =250 W/m 2 ; (b) P 0 (t) 
=250+200cos (0.02t) W/m 2 . 


3.4 Surface Heating 

Heating with a hot plate or pad is a traditional approach to 
retain from pain. Depending upon the temperature and thermal 
properties of heating disk, this approach can be used for cell 
repair or to destroy affected cells. In this section, we analysed 
the thermal behaviour of living tissue subjected to time- 
dependent surface heat flux. Both constant and step heating are 
considered in this study. In constant heating, human tissue is 
heated with a heating pad at a constant rate. In step heating after 
heating for a certain period heat source is removed and allows 
it to cool. Results are calculated at different heat flux and time. 
Here Eq. 3(a) is used for skin surface boundary conditions. 

Constant Heating: The calculated tissue temperature for 
constant surface heating is shown in Fig. 4 (a) at different heat 
flux. And skin temperature at different times, along with the 
distance from the skin surface, is shown in Fig. 4 (b). From Fig. 
4 (a) higher heat flux results in a higher temperature, and 
temperature increases as time increases. At the early stage, 
temperature increases rapidly, but as time increases, increasing 
rate decreases and tends to be stabilized. From Fig. 4 (b) it is 
clear that temperature decreases towards the body core. 



(a) 



(b) 


Fig. 4 Effect of surface heat flux to the skin surface 
temperature response. 
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Step Heating: In this case, after heating for 1200 seconds 
the heat source is removed. In this section, temperature 
distribution at three different locations is obtained over time. It 
is very useful in eye surgery via a single laser pulse due to a 
flash fire, heating using a hot plate for a short period of time 
[10]. The heating power used in this particular type of 
investigation is expressed as 


Q(t)= 


1000 — ,t < 1200 s 

m2 ' 

0 —,t> 1200s 

m2 


(13) 


in surface heating. Also the higher blood perfusion results in 
lower temperature and quick temperature loss (after 1200 
seconds when Q(t)=0). This happens as a higher blood flow rate 
carried away excess heat. Such information is valuable in 
thermal comfort analysis. In practice, the temperature of the 
surrounding fluid temperature and duration should be in the 
safe range. A high temperature or long durable process may 
encounter pain, even burning of the skin. 

3.5 Effect of Heating Frequency and Blood Perfusion 


The transient temperature at three different locations of the 
skin is shown in Fig. 5. Where Q r =0, in Fig. 5. The result also 
carried out for two different value of blood perfusion 
<Db=0.0005 ml/s/ml and CDb=0.004 ml/s/ml. 




(b) 

Fig. 5 Transient temperature at three positions (Qr=0); (a) 
CDb=0.0005 ml/s/ml; (b) cOb=0.004 ml/s/ml. 

Both figures show that as time increases, temperature also 
increases, but after 1200 seconds when the heat source is 
removed, tissue temperature decreases as time passes. 
Moreover, these figures show us the effect of blood perfusion 


The calculated result for different heating frequency and 
blood perfusion is shown in Fig. 6 subjected to sinusoidal 
surface heating. The sinusoidal heating at the skin surface can 
be expressed as 

Q(t)=q 0 +q w cos(coit) (14) 

Where qo and q w are the constant terms, and the oscillation 
amplitude of sinusoidal heat flux and coi represents the heating 
frequency. 



Fig. 6 Effect of heating frequency and blood perfusion on 
sinusoidal surface heat flux. 



Time T(s) 

Fig. 7 Different heating condition and its impact on skin 
surface temperature. 
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In Fig. 6 curves A and B, we use blood perfusion as 
0.0005, where in curve C & D it is 0.004. While in curve A and 
C, we use a heating frequency of value, 0.02 were in B & D, it 
is 0.01. From these figures, we can say that high blood 
perfusion results in lower temperatures where temperature 
response under two different heating frequencies almost 
negligible. 

The calculated tissue temperature result subject to 
simultaneously surface and spatial heating is shown in Curve A 
of Fig. 7. While Curve B represents a single sinusoidal surface 
heating (Q r =0), and curve C represents only sinusoidal spatial 
heating (fi(t)=0). The applied surface and spatial heating are 
Q(t)=1000+500cos(0.02t) W/m 2 and P o (t)=250+200cos(0.02t). 

However, Fig. 8 illustrates the impact of frequencies of 
surface heating. Here the applied surface and spatial heating are 
Q(t)=1000+500cos(0.02t) W/m 2 and P o (t)=250+200cos(0.01t) 
respectively. In curve A, we applied simultaneously sinusoidal 
surface and spatial heating. Curve B represents a single 
sinusoidal surface heating (Q r =0) when curve C represents only 
spatial heating (fi(t)=0). 



In curve A of Fig. 8, the frequency of surface heating and 
spatial heating was the same, and thus, due to the same 
frequency, the resultant temperature appears having the same 
frequency as external heating. However, in curve A of Fig. 8, 
different heating frequency was applied to spatial and surface 
heating; that's why irregular frequency has appeared in tissue 
temperature. 

3.6 Impact Forced Convection Boundary Condition 

In this section, we concentrate on temperature profiles 
under different kinds of surrounding medium classified by their 
temperature. Fig. 9 depicts the tissue temperature distribution 
under different cooling medium temperature. Here force 
cooling significantly reduces the skin surface temperature. 
Moreover, lower cooling medium temperature results in lower 
skin surface temperature. However, the effect of forced cooling 
is negligible for the deep tissues, as shown in Fig. 9, the 


temperature of the cells over approximately x=12 mm line 
remains changeless for different cooling temperature. 



Fig. 9 Influence of cooling medium temperature on tissue 
temperature. 

Using a cooling medium on the skin surface may be a good 
approach during hyperthermia treatment as it can reduce the 
skin temperature even below the core temperature, which may 
result in hypothermia. Hence, concentration should be given 
selecting proper cooling medium. Here the force convection 
coefficient of the cooling medium is considered as 100 W/(m 2 . 
°C). 

3.7 Point Heating 

Treating deep tumors- located at kidney, lung, or rectum- 
it is very difficult to adopt surgical treatment. In such a case, 
point heating can be an alternative to surgery due to its ability 
to treat a tumor with a defined volume. In such a case, the total 
heating power is deposited at the tumor site with the help of a 
microwave probe, radio-frequency probe etc. In this heating 
type, the target region is heated more than 50°C within few 
minutes. This heating is very beneficial in the case of thermal 
ablation when a target tissue is destroyed, injecting thermal 
energy at the tumor site [5],[6]. There is an inverse relationship 
between elevated temperature and exposure duration. For the 
same amount of tissue necrosis, the high temperature needs low 
exposure duration. On the contrary, the low temperature needs 
high exposure duration. So for effective treatment, we need to 
know the required temperature and exposure duration 
precisely. Moreover, in some cases, to protect the skin surface 
cells from excess heat, the cooling medium is used on the skin 
surface during treatment, which is a very efficient approach to 
reduce the skin surface temperature. In this investigation, we 
calculate the tissue temperature at different times, cooling 
medium properties, and heating power. To deposit total heating 
power at the desired site, we use the expression of the external 
heating as [12],[13]. 

Q r (x,t)=P i (t)8(x-x 0 ) (15) 


66 














M. Sannyal and A.M.M. Mukaddes /JEA VoL 01(02) 2020, pp 61-69 


Where Pi(t) is the point, heating power, and 8(x-xo) is the 
Dirac delta function. It has a value 1 at our desired point (xo), 
and at all other points, its value is 0. That’s why it allows 
depositing total heating power at the tumor site. Where xo is the 
distance of tumor site from the skin surface. Here we consider 
the distance of the tumor site from the skin surface is 21 mm 
(xo). In this case, convection boundary condition is applied 
(hf=100 W/(m 2 .°C) and Tf=15 °C) at the skin surface. 



Fig. 10 Impact of point heating on tissue temperature 
distribution. 


In Fig. 10, temperature distribution at different times is 
shown where point heating with a point heat source of 
Pi(t)=2500 W/m 2 is applied. This figure demonstrates that due 
to point heat source, the position of the maximum temperature 
remains constant at the site of the point source at different 
times. 



Fig. 11 Influence of cooling medium temperature to steady 
state temperature distribution. 


In Fig. 11 the temperature response under different 
temperatures of cooling fluid is analysed where in Fig. 12 
influence of tissue temperature under different heating power 
is shown. Both results are computed for the steady state 
condition. In Fig. 11, point source of Pi(t)=2500 W/m 2 is used. 
Fig. 11 shows that the magnitude and position of the highest 
steady state temperature are changeless at different cooling 
medium temperature. It reduces skin surface temperature 
considerably. 

From Fig. 12 it is clear that higher power of the point heat 
leads to a higher temperature. Moreover, tissue temperature 
sensitivity due to point heating power decreases along with the 
distance from point heat source. 



Fig. 12 Impact of point heating power on steady state 
temperature distribution. 

3.8 Tissue Temperature Fluctuation under Stochastic Cooling 
Medium Temperature 

Earlier, we consider the surrounding fluid temperature as 
constant. However, practically surrounding fluid temperature 
does not remain constant; rather, it fluctuates over time. So it is 
necessary to know the impact of such stochastic behavior. For 
this purpose, we use the following expression for flowing 
medium temperature. 

T e =T^s(t) (16) 

Where s(t) the stochastic variance in T e and Tf is the 
equilibrium value if the environmental temperature. This 
variance gives the environmental temperature a stochastic 
value. We assume 

8 t =M0.05-Oi) (17) 

Where 8 and t is the stochastic variance in environmental 
temperature and the discrete-time, respectively; a is the 
random number between 0 and 1. 

Fig. 13 and Fig. 14 depict the influences of variance in 
environmental temperature with different convection 
coefficient between the cooling medium and skin surface. Fig. 
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13 (a) and Fig. 14 (a) shows the tissue temperature fluctuation 
over time where Fig. 13 (b) and Fig. 14 (b) shows the 
fluctuation of stochastic variance. These figures demonstrate 
that due to irregular cooling medium temperature, the tissue 
temperature fluctuates within a certain range. Moreover, the 
frequency of the tissue temperature is much smaller than that 
of the stochastic variance. It may be noticed that as convection 
coefficient gets larger the temperature fluctuation magnitude 
also increases slightly. 


0.4 



°- 4 0 25 50 75 100 125 150 

Time t(s) 

(a) 



Fig. 13 Impact of stochastic temperature variance on tissue 
temperature (hf=100 Wm' 2 ). 

3.9 Tissue Temperature Fluctuation Due to Stochastic Heating 

In this section, we will analyze about stochastic heating, 
which may be encountered for biological rhythm or stochastic 
external heating in hyperthermia treatment. This case 
corresponds to the spatial heating of the following type. 

Qr=Q'm(t) (18) 


Where Q' m (t) is the stochastic variance in metabolic rate, and in 
the initial state, the metabolic rate was considered as constant 
QM. Here we assumed that 

Q'm(t)^ q (0.5-ai) (19) 



0 25 50 75 100 125 150 


Time t(s) 

(a) 



0 25 50 75 100 125 150 

Time t(s) 

(b) 


Fig. 14 Impact of stochastic temperature variance on tissue 
temperature (hi—25 Wm' 2 ). 

Where Q' m and t is the stochastic variance in metabolic 
heat generation and the discrete-time respectively; Gi the 
random number between 0 and 1 and t is the discrete stochastic 
variance in environmental temperature. And X q is a constant, 
which was regarded as Q m /10 in this study. 

The calculated results are shown in Fig. 15 which indicates 
that due to stochastic heating the temperature fluctuates within 
a small range ±0.1 . Fig. 15 (a) shows the tissue temperature 
fluctuation over time where Fig. 15 (b) shows the fluctuation 
of stochastic variance. Section 4.8 and Section 4.9 clearly 
indicates that the biological body, tends to keep its temperature 
balance. 
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(a) 
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(b) 

Fig. 15 Influence of variance in metabolic rate on tissue 
temperature. 

4. Summary 

In this paper, a one dimensional Finite Element Model was 
developed to know the temperature profile inside the human 
tissue subject to numerous heating pattern i.e., spatial heating, 
surface heating, point heating, and stochastic heating. Effect of 
heating frequency, blood perfusion, and scattering coefficient 
are also discussed briefly. Which can be used in parameter 
estimation. It is found that for destroying a target cell point, 
heating is more suitable than other heating as it increases the 
temperature of the target region, and it has a relatively low 
impact on nearby unaffected cells. Moreover, a higher 
scattering coefficient leads to a higher temperature, where 
higher blood perfusion leads to lower temperate. As heating 
apparatus such as laser or microwave may have different power 
and scattering coefficient. The results obtained in this paper can 
be used to select a suitable apparatus. During treatment, 
fluctuation of environmental fluid temperature may be out of 
consideration as its impact on tissue temperature is almost 
negligible, as shown in stochastic heating. The different heating 


styles used for investigation in this study are generally carried 
out in clinical trials. Hence results described in this paper could 
be beneficial to predict the treatment outcome before the 
treatment. This will help to detect the possible risk as well as 
increasing the effectiveness of the treatment. Moreover, the 
developed Finite Element Model and computer code can be 
further use to solve more practical bio-heat transfer problems. 
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